{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 1\n",
    "\n",
    "**Due 2/4/2020 on gradescope.**\n",
    "\n",
    "## References\n",
    "\n",
    "+ Lectures 1-5 (inclusive).\n",
    "\n",
    "## Instructions\n",
    "\n",
    "+ Type your name and email in the \"Student details\" section below.\n",
    "+ Develop the code and generate the figures you need to solve the problems using this notebook.\n",
    "+ For the answers that require a mathematical proof or derivation you can either:\n",
    "    \n",
    "    - Type the answer using the built-in latex capabilities. In this case, simply export the notebook as a pdf and upload it on gradescope; or\n",
    "    - You can print the notebook (after you are done with all the code), write your answers by hand, scan, turn your response to a single pdf, and upload on gradescope.\n",
    "\n",
    "+ The total homework points are 100. Please note that the problems are not weighed equally.\n",
    "\n",
    "**Note**: Please match all the pages corresponding to each of the questions when you submit on gradescope. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context('talk')\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Student details\n",
    "\n",
    "+ **First Name:**\n",
    "+ **Last Name:**\n",
    "+ **Email:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1\n",
    "\n",
    "This exercise demonstrates that probability theory is actually an extension of logic. Assume that you know that ``A implies B\". That is, your prior information is:\n",
    "$$\n",
    "I = \\{A \\implies B\\}.\n",
    "$$\n",
    "Please answer the following questions in the space provided:\n",
    "\n",
    "A. (4 points) $p(AB|I) = p(A|I)$.<br>\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B. If $p(A|I) = 1$, then $p(B|I) = 1$.<br>\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C. If $p(B|I) = 0$, then $p(A|I) = 0$.<br>\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "D. B and C show that probability theory is consistent with Aristotelian logic. Now, you will discover how it extends it. Show that if B is true, then A becomes more plausible, i.e.\n",
    "$$\n",
    "p(A|BI) \\geq p(A|I).\n",
    "$$\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "E. Give at least two examples of D that apply to various scientific fields. To get you started, here are two examples:\n",
    "\n",
    "- $A$: It is raining. $B$: There are clouds in the sky. Clearly,  $A \\implies B$. $D$ tells us that if there are clouds in the sky, raining becomes more plausible.\n",
    "    \n",
    "- $A$: General relativity. $B$: Light is deflected in the presence of massive bodies. Here  $A \\implies B$. Observing that $B$ is true makes $A$ more plausible.\n",
    "\n",
    "**Answer:**\n",
    "\n",
    "- $A$: <br><br>\n",
    "  $B$: <br><br>\n",
    "  \n",
    "- $A$: <br><br>\n",
    "  $B$: <br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F. Show that if $A$ is false, then $B$ becomes less plausible, i.e.: \n",
    "$$\n",
    "p(B|\\neg AI) \\leq p(B|I).\n",
    "$$\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "G. Can you think of an example of scientific reasoning that involves F? For example:\n",
    "$A$: It is raining. $B$: There are clouds in the sky. F tells us that if it is not raining, then it is less plausible that there are clouds in the sky.<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "H. Do D and F contradict Karl Popper's [principle of falsification](https://en.wikipedia.org/wiki/Falsifiability), \"A theory in the empirical sciences can never be proven, but it can be falsified, meaning that it can and should be scrutinized by decisive experiments.\"<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2\n",
    "\n",
    "Disclaimer: This example is a modified version of the one found in a 2013 lecture on Bayesian Scientific Computing taught by Prof. Nicholas Zabaras.\n",
    "I am not sure where the original problem is coming from.\n",
    "\n",
    "We are tasked with assessing the usefulness of a tuberculosis test.\n",
    "The prior information I is:\n",
    "\n",
    "> The percentage of the population infected by tuberculosis is 0.4\\%.\n",
    "We have run several experiments and determined that:\n",
    "+ If a tested patient has the disease, then 80\\% of the time the test comes out positive.\n",
    "+ If a tested patient does not have the disease, then 90\\% of the time the test comes out negative.\n",
    "\n",
    "To facilitate your analysis, consider the following logical sentences concerning a patient:\n",
    "\n",
    "> A: The patient is tested and the test is positive.\n",
    "\n",
    "> B: The patient has tuberculosis.\n",
    "\n",
    "A. Find the probability that the patient has tuberculosis (before looking at the result of the test), i.e., $p(B|I)$. This is known as the base rate or the prior probability.<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B. Find the probability that the test is positive given that the patient has tuberculosis, i.e., $p(A|B,I)$.<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C. Find the probability that the test is positive given that the patient does not have tuberculosis, i.e., $p(A|\\neg B, I)$.<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "D. Find the probability that a patient that tested positive has tuberculosis, i.e., $p(B|A,I)$.<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "E. Find the probability that a patient that tested negative has tuberculosis, i.e., $p(B|\\neg A, I)$. Does the test change our prior state of knowledge about about the patient? Is the test useful?<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F. What would a good test look like? Find values for \n",
    "$$\n",
    "p(A|B,I)= p(\\text{test is positive} |\\text{has tuberculosis},I),\n",
    "$$\n",
    "and\n",
    "$$\n",
    "p(A| \\neg B, I) = p(\\text{test is positive}|\\text{does not have tuberculosis}, I),\n",
    "$$\n",
    "so that\n",
    "$$\n",
    "p(B|A, I) = p(\\text{has tuberculosis}|\\text{test is positive}, I) = 0.99.\n",
    "$$\n",
    "There are more than one solutions. How would you pick a good one? Thinking in this way can help you set goals if you work in R\\&D. If you have time, try to figure out whether or not there exists such an accurate test for tuberculosis<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3\n",
    "Let $A$ and $B$ be independent conditional on $I$. Prove that: \n",
    "$$\n",
    "A \\perp B | I \\iff p(AB|I) = p(A|I)p(B|I).\n",
    "$$\n",
    "Hint: Use the fact that $A \\perp B|I$ means that $p(A|B, I) = p(A|I)$ and $p(B|A, I) = p(B|I)$.<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4\n",
    "Let $X$ be a continuous random variable and $F(x) = p(X \\leq x)$ be it's cumulative distribution function. Using only the basic rules of probability, prove that:\n",
    "\n",
    "A. The CDF starts at 0 and goes up to 1: \n",
    "$$\n",
    "F(- \\infty) = 0 \\ \\text{ and } \\ F(\\infty) = 1.\n",
    "$$\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B. $F(x)$ is a monotonically increasing function of $x$, i.e., \n",
    "$$\n",
    "x_1 \\leq x_2 \\implies F(x_1) \\leq F(x_2).\n",
    "$$\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C. The probability of $X$ being in the interval $[x_1, x_2]$ is:\n",
    "$$\n",
    "p(x_1 \\le X \\le x_2|I) = F(x_2) - F(x_1).\n",
    "$$\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 5\n",
    "Let $X$ be a random variable. Prove that:\n",
    "$$\n",
    "\\mathbb{V}[X] = \\mathbb{E}[X^2] - \\left(\\mathbb{E}[X]\\right)^2.\n",
    "$$\n",
    "**Proof:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 6\n",
    "\n",
    "Hint: Before attempting this example, make sure you understand the Lecture 5 examples.\n",
    "You basically have to repeat the same procedure.\n",
    "\n",
    "The [San Andreas fault](https://en.wikipedia.org/wiki/San_Andreas_Fault) extends through California forming the boundary between the Pacific and the North American tectonic plates.\n",
    "It has caused some of the major earthquakes on Earth.\n",
    "We are going to focus on Southern California and we would like to assess the probability of a major earthquake, defined as an earthquake of magnitude 6.5 or greater, during the next ten years."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A. The first thing we are going to do is go over a [database of past earthquakes](https://scedc.caltech.edu/significant/chron-index.html) that have occured in Southern California and collect the relevant data. We are going to start at 1900 because data before that time may are unreliable.\n",
    "Go over each decade and count the occurence of a major earthquake (i.e., count the number of organge and red colors in each decade). We have done this for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-4976fee020d5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m eq_data = np.array([\n\u001b[0m\u001b[1;32m      2\u001b[0m     \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;31m# 1900-1909\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;31m# 1910-1919\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;31m# 1920-1929\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;31m# 1930-1939\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "eq_data = np.array([\n",
    "    0, # 1900-1909\n",
    "    1, # 1910-1919\n",
    "    2, # 1920-1929\n",
    "    0, # 1930-1939\n",
    "    3, # 1940-1949\n",
    "    2, # 1950-1959\n",
    "    1, # 1960-1969\n",
    "    2, # 1970-1979\n",
    "    1, # 1980-1989\n",
    "    4, # 1990-1999\n",
    "    0, # 2000-2009\n",
    "    2 # 2010-2019 \n",
    "])\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar(np.linspace(1900, 2019, eq_data.shape[0]), eq_data, width=10)\n",
    "ax.set_xlabel('Decade')\n",
    "ax.set_ylabel('# of major earthquakes in Southern CA');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B. The right way to model the number of earthquakes $X_n$ in a decade $n$ is using a Poisson distribution with unknown rate parameter $\\lambda$, i.e.,\n",
    "$$\n",
    "X_n | \\lambda \\sim \\operatorname{Poisson}(\\lambda).\n",
    "$$\n",
    "Here we have $N = 12$ observations, say $x_{1:N} = (x_1,\\dots,x_N)$ (stored in ``eq_data`` above).\n",
    "Find the *joint probability* (otherwise known as the likelihood) $p(x_{1:N}|\\lambda)$ of these random variables.<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C. The rate parameter $\\lambda$ (number of major earthquakes in per ten years) is positive. What prior distribution should we assign to it.\n",
    "A suitable choice here is to pick a [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution), see also [the scipy.stats page for the Gamma](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html).\n",
    "We write:\n",
    "$$\n",
    "\\lambda \\sim \\operatorname{Gamma}(\\alpha, \\beta),\n",
    "$$\n",
    "where $\\alpha$ and $\\beta$ are positive *hyper-parameters* that we have to set to represent our prior state of knowledge.\n",
    "The PDF is:\n",
    "$$\n",
    "p(\\lambda) = \\frac{\\beta^\\alpha \\lambda^{\\alpha-1}e^{-\\beta \\lambda}}{\\Gamma(\\alpha)},\n",
    "$$\n",
    "where we are not conditioning on $\\alpha$ and $\\beta$ because they should be fixed numbers.\n",
    "Use the code below to pick some some reasonable values for $\\alpha$ and $\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEdCAYAAAA1s6EDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3wb15Xo8R/A3kVRonqXfSjbklXdey+JI3c5WTt1kzjeFO9L74mzcZxN27y8JN44bokd9xYn7k2ukWQ1qx31TnWxiL3g/XEHEgSRIkASGBA8388HnyGm4eACxJk7c+feQCgUwhhjjEmUoN8BGGOMSW+WaIwxxiSUJRpjjDEJZYnGGGNMQlmiMcYYk1CZfgeQglpxCbjG70CMMaaPKAba6SSnBKx58xHaQ6FQoDvFEgi4aX8vUisHx8rBsXJw0rkcAgEIBAIhOjlLZjWaI9WEQpTs3Xsg7g1LSvIAqK5u6O2Y+hQrB8fKwbFycNK5HMrKCgkEOj8LZNdojDHGJJQlGmOMMQmVsqfORGQqMB8Yp6pbj7JeIXAHcDVQCMwFvqyqa5ISqDHGmKNKyRqNiAjwLLElwoeBa4FvADcBI4DXRKQkcREaY4yJVUrVaEQkE/gs8DOgJYb1zwAuAy5V1ee9eW8CG4DP42o6xhhjfJRqNZozgJ8Dv8TVULpyEVALvBSeoaq7gTdwCcgYY4zPUi3RrATGq+qPcDdOdqUCWKuqbVHz1wLS28F1ZdXGfeypSr+mi8YY0xMpdepMVXfGuUkJHd/BX4u7U7VbAoFDbd5jtXZLFd+5812GDSrgt7eeTTAY6O7L93mZmRlA/GWYbqwcHCsHJ53LIdDFz12q1WjiFQA6us82gOsOIWkK8rIAqNxTx6pN+5P50sYYk9JSqkbTDdXA+A7mF3nLuiUUiv/u3fysIOOGF7Nhew2vL9jCiIHpd9QSq3S+AzoeVg6OlYOTzuXg9QzQqb5eo1FgvIhEv8WJ3rKkOn3KcAAW6C7a29OwQyNjjOmGvp5oXgQGABeEZ4jIYOAs4OVkB3Pa5GEAVNc1o1uqkv3yxhiTkvpUohGRwSJyiogUA6jqXOB14CER+bSIXIlLMFXAH5Id35CB+Uwc6e4Tnb9qV7Jf3hhjUlKfSjTA5cC7wPSIeVcBzwC/AO4FtgLnq6ovV+TDp8/e1120tSe1PYIxxqQkG4/mSFXt7aFuDxOwp6qBz93xKgD/5/qpHD9uYG/Hl/LS+aJnPKwcHCsHJ53LoayskGAwUI27lHGEvlajSXmDBuQxcYQ7fTZvZby3BRljTPqxRJMAsyrKAVi4ejetbXb6zBjTv1miSYBZk8oJAHWNrSxbv8/vcIwxxleWaBJgQGEOFWNKAXhvxQ6fozHGGH9ZokmQU44fAsDiNXtoaIqlf1BjjElPlmgSZMax5WRmBGlubWfh6t1+h2OMMb6xRJMg+bmZTJ1YBsB7K6z1mTGm/7JEk0AnHzcUgBUb91F9oMnnaIwxxh+WaBJoyoQy8nMyCYVg3krrksYY0z9ZokmgrMwgMysGA9b6zBjTf1miSbBTvNNnGyprqdxb53M0xhiTfJZoEuzY0QMoK84B4J1lVqsxxvQ/lmgSLBgIcOoJbpyad5btsAHRjDH9jiWaJDj9BHf6bH9tEys3+zJ6gTHG+MYSTRIMGZh/sEfndz6o9DkaY4xJLks0SXLaZFereX/1buuSxhjTr1iiSZKTKsrJygzS3NLOArV7aowx/YclmiTJz81i2jGDAHjnA2t9ZozpPyzRJNHpk13rM91Sxa6q9BvO1RhjOmKJJomOHzuQ0iJ3T81bS61RgDGmf7BEk0TBYOBgrebtDyppa7dhno0x6c8STZKdOcUlmv21TTbMszGmX7BEk2SDB+QxyRvmee6S7T5HY4wxiWeJxgdnnTgcgKXr9to4NcaYtGeJxgfTjx1EQW4mbe0h62jTGJP2LNH4ICszg1O9/s/mLtlOKGQdbRpj0pclGp+cNcWdPtu5vwHdXOVzNMYYkziZfgcQTURuAL4LjAc2Arer6v1HWX8w8HPgYiAXeAe4VVXXJD7a7htZXsiE4cWs217D64u3UeE1EDDGmHSTUjUaEbkWeAB4EZgNvA7cJyLXdLJ+AHgSuBT4JnAjMBR4TURS/pf7nGkjAHhfd1Nd1+xzNMYYkxgplWiA24FHVPVWVX1BVW8GHgFu62T9Y4DTga+r6v2q+g/gOmAEcEVSIu6BWRXlBxsFvLXUmjobY9JTyiQaERkPTAAej1r0GFAhIuM62CzXm9ZGzAvfBVnWuxH2vuysjIM9BbyxeLuNvmmMSUspk2iACm+qUfPXelOJ3kBVlwKvAd8XkQrves1vgQPAU4kKtDedPdU1CthT3ciyDdZTgDEm/aRSY4ASb1oTNT9cWynuZLubgReAld7zJmC2qq7vbiCBAJSU5MW9XWZmBhDftiUleUyeUMYH6/by1rJKzpw+Mu7XTTXdKYd0ZOXgWDk46VwOgcDRl6dSjSYcavT5o/D8I3qgFJFJwHvAbuBKXMuzvwOPi8iZCYqz11188hgAFq7axW4bPsAYk2ZSqUZT7U2jay5FUcsj3epNL1LV/QAi8hLwJvBrYGZ3AgmFoLo6/h/88JFKvNseO6KYksJsqg808/e567j67Alxv3Yq6W45pBsrB8fKwUnncigrKzxqrSaVajThazMTo+ZPjFoeaQywIpxkAFQ1BLwFHN/rESZIZkaQc6a6ps5vLN5Oc0ubzxEZY0zvSZlEo6prgQ1A9D0zVwNrVHVzR5sBJ3Rwz8wpuJs9+4xzpg4nIxjgQEML/1q50+9wjDGm16TSqTOAHwP3iMh+4FncvTDXAXPgYC8AE3C1mBrgV8C/AS+IyM+AeuAm4OzwNn1FSWEOsyaV897ynbyyYCtnTB5GoKsrbMYY0wekTI0GQFXvBT6Pu6j/FHAOcJOqPuytcjnwLjDdW38j7obNHcC9wEPAKODCiG36jPNnuBZnm3cdYM3Wji5JGWNM3xOwnoOPUNXeHirZu/dA3Bv2xsW+2+6bz4bKWmZWlPOF2Sd0ez9+SueLnvGwcnCsHJx0LoeyskKCwUA1MKCj5SlVozFwwYxRACzU3eyrafQ5GmOM6TlLNClmZkU5xQXZtIdCvLJwq9/hGGNMj1miSTFZmUHO83p1fmPRdhqbW32OyBhjesYSTQo6Z/oIsjKD1De18tbSSr/DMcaYHrFEk4KK87M53Rvq+aUFW6xXZ2NMn2aJJkVdOMs1Cthd1ciiNbt9jsYYY7rPEk2KGlZWwIkT3JA6L8zb4nM0xhjTfZZoUtjFJ40GYO22atZtsxs4jTF9kyWaFCajBzBmiOu8+rl/ddTVmzHGpD5LNCksEAhw6SmuVrNw9W6276nzOSJjjImfJZoUN1PKKR/guq543mo1xpg+yBJNigsGA1zi1WreXb7DuqUxxvQ5lmj6gNNPGEpJQTZt7SFenG8t0IwxfYslmj4gKzODi7z7at5YvJ0DDS0+R2SMMbGzRNNHnDNtBHk5mTS1tPHyAqvVGGP6Dks0fUReTibnTXedbb68YCv1jdbZpjGmb7BE04dcNGsUOVkZ1De12hACxpg+wxJNH1KUn32wVvPivM00NFmtxhiT+izR9DEXnzSa7MwgdY2tvLZom9/hGGNMlyzR9DHFBdmc4w2M9vy/NtPU3OZzRMYYc3SWaPqgS04eTWZGkAMNLVarMcakPEs0fdCAwhzOnjocgH++t8mGezbGpDRLNH3U5aeOITvT1WpeXmAt0IwxqcsSTR81oDCH82aMBNy1mvpG6y3AGJOaLNH0YZeePJqcbHdfjfWBZoxJVZZo+rCi/GwunOn6QHtx/hZq65t9jsgYY45kiaaPu+SkUeTnZNLY3GajcBpjUlKm3wFEE5EbgO8C44GNwO2qev9R1g8C3wI+DQwD1gL/paoPJT5a/+XnZnHJyaN5Yu56Xl6wlQtmjGRgca7fYRljzEEpVaMRkWuBB4AXgdnA68B9InLNUTb7DfA94HfAh4D3gAdF5NLERps6Lpw5ipKCbFrb2nnqrQ1+h2OMMYdJtRrN7cAjqnqr9/wFERkI3AY8Fr2yiEwAbgE+q6p/9ma/IiLHApcAzyUhZt/lZGfwkTPGcf8LytsfVHLxrFGMGFzod1jGGAOkUI1GRMYDE4DHoxY9BlSIyLgONpsN1AOHnVpT1bNV9csJCTRFnTFlGEMG5hMKweNvrPc7HGOMOShlEg1Q4U01av5abyodbDPFW/9CEVkiIq0iskZErk9UkKkqMyPI1WeNB2Dx2j2s3lLlc0TGGOOk0qmzEm9aEzW/1psWd7DNYGA0cDfuOs0G4DPAQyKyS1Vf604ggQCUlOTFvV1mZgbQvW17w3knjealBVtZs7WKx+eu5/abTyMQCCQ9Dr/LIVVYOThWDk46l0NXPzOpVKMJhxrqZH57B9tk45LNZ1T1T6r6MnADsAT4YSKCTGWBQIAbL3MVwzVbqnhryXafIzLGmNSq0VR70+iaS1HU8ki1QBuulRoAqhoSkZdwNZtuCYWguroh7u3CRyrd2ba3jByYz4xjB/P+6t3c/9xKKkaWkJ2VkdQYUqEcUoGVg2Pl4KRzOZSVFR61VtOtROO16joeKMfVQHYDy1R1TXf25wlfm5kIfBAxf2LU8khrcLWyLCDytvhsjqwZ9RvXnjuBxWv3sK+miRfnb+FDp431OyRjTD8W86kzEZkkIv8jItuBlbjWYH8A/uj9vUpEtovIb0RkUryBqOpa3DWW6HtmrgbWqGpHt70/jzu1dl1EnJm4ps1vxhtDuigvzeeCma7DzX+8t4nqA00+R2SM6c+6rNF496rcAVwJNOB+wN8F1gF7cT/0A3E1j1Nwp6y+KCJPAN9Q1Xja2v4YuEdE9gPPAlfgksgcL5bBuCbQK1S1RlVfFZF/Ar8VkUJgNfAFYBzw0TheN+18+LSxvP3BDg40tPD43PV86rK4c78xxvSKWE6drcCdyvoE8ISq1h1tZREpwNVKvuRtG3N/KKp6r4jkAF/FJaz1wE2q+rC3yuXAPcC5uF4D8F7rx8A3cQlvEXChqr4f6+umo/zcLGafOY6/vriat5dWcu60EYwb1lHDPWOMSaxAKHT0Sxki8hFVfbo7O+/Jtj6qam8PlezdeyDuDVPtYl9bezs/umcBW3cfYNywYr5z0wyCSWjunGrl4BcrB8fKwUnncigrKyQYDFQDAzpa3uU1mp4kij6YZNJKRjDIxy48BoANlTW8vbTS54iMMf1RKt1HYxJARpdyynFDAHjsjXU2EqcxJuks0fQD1547kZzsDGrrW3jqTevd2RiTXL1yw6bXpLgCmOw9TlDVK3pj36bnSotyuOK0sTz6+jpeWbiV0ycPY8zQoq43NMaYXhB3ohGRMRxKKJOBE3AdXob3FQDs/EyKuXDWKN5etoPte+q4/4VVfOfGmQSDye8HzRjT/8RyH81oXHPjGbikEjnQSQDYBbwBLI14rOj1SE2PZGYEueli4WcPLGRDZS2vLdrG+TNG+h2WMaYfiKVG8zBwMnAAWIbr3mU6UAlcrarvJS4805uOHTWAM6YM462llTwxdx3Tjx1MaVGO32EZY9JcLI0BpgF/AkpV9TRVnQn8B1AAvCwi/ykidg6mj7ju3IkU5mXR0NTGQ6/0pGs6Y4yJTSyJ5n3geVVtC89Q1d8Dk3C9Jv8CeE9EJicmRNObCvOyuO5c10/p/FW7WLRmt88RGWPSXSw3bJ6uqk92ML9SVa8CrgJGAAtE5Ccikp2AOE0vOn3yUCaNKQXgLy8o9Y2tPkdkjElnPb6PRlWfwtVu7sL1N7ZERM7o6X5N4gQCAT5xaQXZWUGqDjTz6Otru97IGGO6qVdu2FTVWlW9BTgDaOVQh5cmRQ0ekMfVZ00A4I3F21m5ab/PERlj0lWv9gzgtUCbBny/N/drEuP8GSOZMML16HzvcytpbLZTaMaY3tdlohGR8+PZoaq2qupPvW0v6G5gJvGCwQCfuHQSmRkBdlc18ujr6/wOyRiThmKp0TwvIq+KyIdEpMvB50UkS0SuFJE3gH/2PESTSCMGFXDlWeMBeG3hNpZv2OdzRMaYdBPLDZvTgF8BzwB7ROQlYB5uhM19HBph8xjcCJvnec9fAKYmIGbTyy6eNZpFq/ewdls1d/9zJbd9+iTyc7P8DssYkya6HPgsTEROxQ2T/BFcNzTRGwaAGuAJ4A+qOr8X40ymtBn4LB4799fzg7vn0dzSzmknDOUzHzqu2/vqy+XQm6wcHCsHJ53LoauBz2LuVFNV3wXe9U6fzQCOAwbjEs5uXPc0i1S1vcdRm6QbUprPtedM5IGXVvPOsh1MnTiImRXlfodljEkDcffe7PUQMM97mDRy7vQRLF67h+Ub9nHf86sYP7yYgcW5fodljOnjujUejYgMAa4DxuI621wIvKSq9b0Xmkm2YCDApy+fxPf/PI8DDS3c9ewKvjpnmg0nYIzpkbjvoxGRM4G1wG+AW4HvAU8Cm0TkS70bnkm2AYU5fPKyCgBWba7ihXmbfY7IGNPXdeeGzV94008Bo3G1mjm4Vmi/EZEHeyc045dpxwzm3GkjAHhi7nrWb6/xOSJjTF/WnURzAvBrVb1PVbeq6mZVfURVTwH+HbheRP6jd8M0yXbdeRMZPqiAtvYQf3x6GXWNNmiqMaZ7upNoaoEOz6eo6p+Bh4DP9yQo47+crAxu/sjxZGcG2VPdyN3/WEmsTeGNMSZSdxLNa8BlXSyf0L1wTCoZMbiQj114LACL1uzh5fe3+hyRMaYv6k6i+V/gNBH5cifLxwLbux2RSSlnTBnGqccPBeCRV9eybnu1zxEZY/qa7jRvfgU3FMCvRORK3Dg073vzzga+DHyj1yI0vgoEAtx48bFs3FFD5d56fv/kMn7wiVkUF9j4dsaY2HSnRnMbrrPMLcBZwP24XgFWAXcCClSJyAkiEnciE5EbRGS5iDSIyEoRuSmObUeJSLWIfDfe1zWdy83O5JYrJ5OTncH+2ib++PQy2tqtAwhjTGziTjSq+gNVna2qY3GdZ54PfBX4K7AcmOL9vQQ4ICJLROQvsexbRK4FHgBeBGbjBlC7T0SuiWHbAHA3UBzvezJdGz6ogE9dNglw99c8MXe9zxEZY/qKbvUMEKaqVbiL/6+F54lINq4J9FRcz8/TgA/HuMvbgUdU9Vbv+QsiMhBXi3qsi21vBipij97Ea1ZFOetPGsUL87bw3HubGTu0mFnWH5oxpgs9SjQdUdVmXJc0C+PZTkTG41qrfStq0WPAdSIyTlU3HGXbO4BrgefiDtrE7JpzJrBpRy2rNlfx52dXUD4gjzFDi/wOyxiTwnp1KOceCtdGNGr+Wm8qHW0kIkHgXlxN6PnEhGbCMoJBbp59AoNKcmlubef/PrGU6rpmv8MyxqSwXq/R9ECJN43u76TWm3Z27eUrwHhiPz3XpUDg0NgR8cjMdAOQdmfbvqSkJI9vfXwW3/njO+yraeLOZ5bzw8+cTJb3/vtLOXTFysGxcnDSuRwCXfS7m0o1mnCoHQ2oBnBEMycREeAnwL+rqt3gkURjhxXzpevcAKqrNu3nj08us54DjDEdSqUaTThRRNdciqKWA+ANwHYf8CjwUlRT6qCIZKpqa3cCCYW6NwpeOo+g15GKkSVcedZ4npy7ntcXbqW0IIsPnz6u35VDZ6wcHCsHJ53Loays8Ki1mlSq0YSvzUyMmj8xannYKOBk4CagJeIB8KOIv00CfejUMZw+2fUc8OSbG3hvxQ6fIzLGpJqUSTSquhbYAETfM3M1sEZVozvy3A7M6uAB8IeIv00CBQIBPn5JBRWj3VDhd/9jJSs27PM5KmNMKkmlU2cAPwbuEZH9wLPAFbiRPOcAiMhgXBPoFapaAyyI3oG7bMN2VT1imUmMzIwgt1w1mZ/+5X0q99bzs/vn85PPnUZJXqp9vYwxfkiZGg2Aqt6LG2LgYuAp4BzgJlV92FvlcuBdYLof8ZnOFeRmceu1J1JSkE1dYyu33TOPfTWNfodljEkBAWspdISq9vZQyd69B+LeMJ0v9sVq885afv7gIuqbWhlWls+3/m0GhXlZfoflC/s+OFYOTjqXQ1lZIcFgoBoY0NHylKrRmL5v9JAivn7jDDIzglTurec3jy6hoalbjf+MMWnCEo3pdZMnDOLL108lEID122v43RMf0NLa5ndYxhifWKIxCXHa5GF8/BLXq9DKTfv5w1PLaW2zoQWM6Y8s0ZiEOevE4cw5/xgAFq/dw5/+vsLGsTGmH7JEYxLqolmjmH3GOADmr9rFXc+utGRjTD9jicYk3IdPH8vlp44B4F8rdvLnf6ykvd1aOxrTX1iiMQkXCAS46qzxXHrKaADeW76Tu/5hp9GM6S8s0ZikCAQCXHP2BC456VCyufOZFdZAwJh+wBKNSZpAIMC15044eBptwapd/P7JZbS0WrIxJp1ZojFJFT6NNvtM10Bg8do9/PaxJTQ2202dxqQrSzQm6QKBAFecPo5rz50AwPKN+/nFQ4s50GAjOxiTjizRGN9cevIYbrpECOB6EPjZAwvZX9vkd1jGmF5micb46pypI7h59glkBANs31PHf/1lAdt2x9+hqTEmdVmiMb6bWVHOV649kZzsDPbVNHH7Xxeim/f7HZYxppdYojEp4fhxA/nmR6dTUpBNfVMrv3x4sQ0LbUyasERjUsaYoUV856YZDCvLp7UtxP8+s4Kn39qAjZlkTN9micaklEEleXz7xhlUjHbjJz391gbufGY5zS02zIAxfZUlGpNyCnKz+M/rp3LWicMAmLdyF3c8uMhapBnTR1miMSkpMyPIxy+p4PrzJhIANlTW8KN757Nma5XfoRlj4mSJxqSsQCDAxSeN5svXnkh+TiY1dc38/MFFvLZom123MaYPsURjUt6UCWV87xMzGTGogLb2EH95Qbnr2ZU0Ndt1G2P6Aks0pk8YUprPt2+cwUwZDMC7y3fwk/sXULm3zufIjDFdsURj+oy8nExunn0Cc84/hoxggG176vjxvQt4+4NKv0MzxhyFJRrTpwQCAS6aNYpvfHQ6pUU5NLW08ed/rORPf19OQ5P1AG1MKrJEY/qkiSNL+OEnZzF14iAA3l2+kx/dO59126t9jswYE80SjemzivKz+eLVk7nhgmPIzAiwa38Dt/9lIU+9ud6GiTYmhViiMX1aIBDgwpmj+O5NrlVaeyjEM29v5Kd/WWgNBYxJEYFUux9BRG4AvguMBzYCt6vq/UdZfyhwG3ARMBBQ4A5VfbSbIVS1t4dK9u6Nv6v6kpI8AKqrG7r50unBr3JoaW3j8TfW8+L8LYC76fPKs8Zx8azRBIOBpMYC9n0Is3Jw0rkcysoKCQYD1cCAjpanVI1GRK4FHgBeBGYDrwP3icg1nayfAzwPXAh8H7gKeB94xEtYph/JysxgzvnH8LUbpjGoJJfWtnYefW0dt//1fRvjxhgfpVSNRkTWAgtUdU7EvIeBKao6qYP1ZwNPAiep6vyI+c8Bw1R1ajfCsBpND6VCOTQ0tfLY6+t4bdE2ADKCAS47ZQwfOm0MWZkZSYkhFcohFVg5OOlcDn2mRiMi44EJwONRix4DKkRkXAeb1QD/CyyImr/K25fpp/JyMrnxYuFrc6ZSPiCPtvYQf39nIz+4ez4rN9mgasYkU6bfAUSo8KYaNX+tNxVgQ+QCVX0VeDVynohkAZcDy7sbSCBw6OgjHpnekXJ3tk0nqVQOp5w4gmnHDeWxV9fw9Nz17NhXz3//bRFnnjicj182idLi3IS9diqVg5+sHJx0LodAF5dAU6ZGA5R405qo+bXetDjG/dwBHAPc3htBmb4vJyuDj11cwc//4wxkTCkAby7Zzpd+9QbPvLmellZrCm1MIqVSjSacE6MvGoXnH/XXQEQCuCRzK/Dfqvp0dwMJhbp3HjWdz8HGI1XLoTQ/i6/NmcrbSyt59PV1HGho4b5/ruSF9zYx5/xjmDKhrFdfL1XLIdmsHJx0LoeyssKj1mpSqUYTvqU7uuZSFLX8CF7rsweBr+GSzNd7PzyTDoKBAGeeOJyffvYUzp8+kmAgwI599fzm0SX86uHFbN1lrdOM6W2plGjC12YmRs2fGLX8MCJSDLwEXAd8xZKMiUVhXhYfu+hYfvipWUzyTqct27CPH9wzj7v/udJG8zSmF6VMolHVtbiL/dH3zFwNrFHVzdHbiEgG8DRwCjBHVf8n4YGatDJycCFfnTOVL109hWFl+YRC8NbSSr5557s88upaDjS0+B2iMX1eKl2jAfgxcI+I7AeeBa7A1VTmAIjIYFyz5RWqWgN8HjgHuBPYIiKnROwrpKr/SmLspo8KBAJMPWYQkycMZO6SSp5+awM1dc08P28zbyzZxkWzRnPhzFHk56bav4sxfUNK/eeo6r3e9ZavAp8B1gM3qerD3iqXA/cA5+J6Dbjam/857xGpjRR7fya1ZQSDnDttBKcdP5SX39/CP9/bTENTK0+/tYGX5m/hopNGccEMSzjGxCulegZIEdYzQA+lSzkcaGjhhXmbeXnBVppa3LDReTmZnD9jJBfOHElRfvZRt0+XcugpKwcnncuhq54BLNEcyRJND6VbOdTWu9Nor76/7WDCycnK4Oypw7lo1igGdnLTZ7qVQ3dZOTjpXA6WaOJniaaH0rUcauubeXnBVl55fyv13mieGcEAJ00awiUnj2ZUeeFh66drOcTLysFJ53KwRBM/SzQ9lO7l0NDUymuLtvHSgi1UH2g+OH/SmFIunDmKKRPLCAYCaV8OsbJycNK5HLpKNHZV05g45eVkctkpY7hw5ijeW76D5+dtpnJvPSs37Wflpv2UD8jjnGkjuOyMcV1exzGmP7AazZGsRtND/a0c2kMhVmzYx4vzt7Bsw76D87Mzg5xx4nBOPX4I44cVE+iq58E01d++D51J53KwU2fxs0TTQ/25HLbvqeO1hdt4e1kljc1tB+ePKi/krBOHc8rxQyjIzfIxwuTrz9+HSOlcDpZo4meJpoesHNx1nMXr9/HSvM1srDzUIXlmRpDpxw7ijCnDOG7MQF+GmE42+z446VwOlmjiZ4mmh6wcnJKSPM8mUOQAABkYSURBVEKhEEt0F68v2sa8VbtoiqjllBRmc+pxQzn1hKGMHFyQtqfW7PvgpHM5WKKJnyWaHrJycKLLobG5lQWrdvPW0u2s3np4Z+TDBxVw8qRyTjpuCENK85MeayLZ98FJ53KwRBM/SzQ9ZOXgHK0cdlU18N7yHbyzbAe79h++fPSQQmZVlDOzojwtko59H5x0LgdLNPGzRNNDVg5OLOUQCoXYuKOWf63YybyVO6mKuC8HYOTgAqYfO5jpxw5mVHlhnzy9Zt8HJ53LwRJN/CzR9JCVgxNvObS3h1iztYoFq3azYPWuw24GBSgrzuHEiYOYOnEQMrqUrMyUGeXjqOz74KRzOViiiZ8lmh6ycnB6Ug7toRDrt9WwcPVuFq7eza6qw/eRk5XBpDGlTJ5QxuRxAxk0IK9XYk4E+z446VwOlmjiZ4mmh6wcnN4qh1AoxLY9dSxZu4fFa/ewflsN0f+1Qwbmc8LYgRw3thQZXZpSQxnY98FJ53KwLmiM6eMCgQAjBxcycnAhl586lpr6Zpav38fS9XtZtn4vdY2t7NxXz8599byycCuBAIwdWsykMaVUjBnAMSMGkJOd4ffbMP2Y1WiOZDWaHrJycJJRDu3tITbtrGXZ+r0s37CPddtraGs//H86Ixhg7NAijhk1gGNHDeCYkSVJ7Z3Avg9OOpeDnTqLnyWaHrJycPwoh6bmNlZvrWLFxn2s2lzF5p21dPQvPnxQARNHFDNhRAkThpcwtCyfYIJatNn3wUnncrBTZ8b0IznZGUweX8bk8WUA1De2sHprNau3VLFmSxUbd9TS1h5i+546tu+pY+6SSsD1SD1+WBFjhxUzznsMKMzuk82pTeqxRGNMGsvPzWKq1yQaoKmljY2VNazZWs3abdWs317DgYYWGppaWb5xP8s37j+4bXFBNmOHFjF6SBFjhhQyekgRg0pyLfmYuFmiMaYfycnKQEa7lmngWrTtqmpg3bZqNlTWsqGyhs07D9Da1k5NXTNL1+1l6bq9B7fPy8lk5OACRpUXMrK8kJGDChkxuIC8HPspMZ2zb4cx/VggEGBIaT5DSvM57YRhALS2tbNtdx2bdtaycUctm3bUsnX3AVpa22loamXN1mrWRPXVVlacw/BBhQwflM/wsgKGDSpgWFl+vxsSwXTMEo0x5jCZGUHGDC1izNAizjrRzWtrb2fHvga27Kxl6+46tu4+wJZdB9hf2wTA3pom9tY08cH6vYftqzg/ixHlRYwYXEBpYbaX1PIoL80jK9OaXPcXlmiMMV3KCAYZMaiAEYMKDptf39jC1t11bPMaF2zfU8f2vXUHu8+pqW+hZuM+Vm7cd9h2AaC0OIfyAXkM9h6DBuS6v0vyKMrPsmtBacQSjTGm2/JzszjWuz8nUn1jK5X76qjcU09VfTPb99SxdWctu/Y30NzaTgjYV9PEvpomVm2uOmK/2ZlBykpy3aP40GNgcQ6lxbmUFub0mb7ejCUaY0wC5OdmMmG4u0cn8v6R9lCI6gPN7NxXz66qBnZXNbBrv5vurmqgrrEVgObWdir31lO5t77T1yjOz6K0KJfSohwGFOUwoDCbAYU53iObksIcivKzEnZ/kImdJRpjTNIEAwFKi3IoLcqhYkzpEcvrG1vZU93A3upG9lQ3srfGe1Q3sq+2iZq6Qz1a19S3UFPfwqadtUd9vaKCLEoKsikuyKYk302L8rMpLsiiON/9XZSfRVF+ll03ShBLNMaYlJGfm8noXHfvTkdaWtvZX9vI/tqmg499tU1UHWiiKjw90HywG55wDSp6yIXO5GRnUJSXRWHEoyA8zc08+Dw/N5OCXDfNz8kkM8NO4x1NyiUaEbkB+C4wHtgI3K6q9x9l/ULgDuBqoBCYC3xZVdckPlpjTDJlZQYpL82n/Cgjj4ZCIeoaW6mqbaK6rpnqOjet8R7Vdc3U1rdQW++mkX3DNTW30dTcxp7qxrjiysnKcEnHSzx5Od7Ue56bncHAAfnk52QSamsnLyeD3Gw33z0yyc4Kpm0DiJRKNCJyLfAA8D/A88Bs4D4RqVfVxzrZ7GFgFvA1oBb4AfCaiByvqtWdbGOMSVOBQOBgbWRkF+uGQiHqm1oPJp4D9S3UNri/6xpaOdDQwoGGFuoaW6hrdM/rG1tobTu8A7mmljaaWtoONvfuXtwcTDo5WRnkZGeQ601zsg49srODh/7OyiAnK0h2ZsTfWRlkZQbJznR/Z2cGycrMIDMj4FsiS6lEA9wOPKKqt3rPXxCRgcBtwBGJRkTOAC4DLlXV5715bwIbgM/jajrGGNOhQCBAQW4WBblZDB3YeS0pUigUormlnbrGFuobWw9O65u8R2MrDd7fDd7zxuZWmlra3bzGVppa2jrYLzQ0tdHQdOSy3hAAsrKCZGV4ySgjSFbm4Y/hgwq47tyJvX4qMGUSjYiMByYA34pa9BhwnYiMU9UNUcsuwtViXgrPUNXdIvIGLgFZojHG9KpAIOBqGdkZDCyOfbvI1ndt7e00NbfR2NxGQ3MbjU2tNLa00djURmNzK80tbQefN7W0HXze1Oz+bmpp96bu0dLaTlNLO61t7Z2+fghobmn3kmRrh+us2Lifs6YMZ2R5YTxF0qWUSTRAhTfVqPlrvangairR26xV1ehDgLXA9d0NJBA49KWIR6bXYqU726YTKwfHysGxcnCSUQ5t7SGaW9oOJqHmlnaaW9u8eS45Nbe6aUtrOy2tbn5Lq1tvyMB8jps4KO5TbF2tnkqJpsSb1kTND7dd7OjYoaSD9cPbxHGsYYwxfV9GMECe1xghlaRSNOGcGD1MU3h+R3XCQAfrh+d3XofsQijUvcGJ0nlgo3hYOThWDo6Vg5PO5VBWVnjUWk0qNf4OtxCLrokURS2P3qajmktRJ+sbY4xJslRKNOFrMxOj5k+MWh69zXgRic6lEztZ3xhjTJKlTKJR1bW4i/3XRC26Glijqps72OxF3BjVF4RniMhg4Czg5QSFaowxJg6pdI0G4MfAPSKyH3gWuAK4DpgDB5PIBGCFqtao6lwReR14SES+DuwDfghUAX9IfvjGGGOipUyNBkBV78XdaHkx8BRwDnCTqj7srXI58C4wPWKzq4BngF8A9wJbgfNVdT/GGGN8FwiFOmq01a+1h0KhQHeKJdzqor8XqZWDY+XgWDk46VwOgQAEAoEQnVReLNEcqRVXWB3dn2OMMeZIxbhbSjq8HGOJxhhjTEKl1DUaY4wx6ccSjTHGmISyRGOMMSahLNEYY4xJKEs0xhhjEsoSjTHGmISyRGOMMSahLNEYY4xJKEs0xhhjEsoSjTHGmISyRGOMMSahUm08mj5LRG4AvguMBzYCt6vq/b4G5SMRmQrMB8ap6la/40kmEQkCnwW+gPs+7ASeBn6gqrV+xpZM3si3X8aVwyhgNXCHqj7oa2A+E5EngCmqGj2acNqyGk0vEJFrgQdwI37OBl4H7hOR6NFC+wUREdzAdf31QObrwO+Af+C+D78EPg486mdQPvgWbpyo+4APAS8BD4jIdb5G5SMR+TfgSr/jSDbrvbkXiMhaYIGqzomY9zDuqGWSf5Ell4hk4o7kfwa0AAOBUf2pRuMdxe8F/qaqt0TMvx54CJimqov9ii9ZRCQLV5N7QFW/GDH/dSBDVc/0Kza/iMhwYBlQBzRZjcbETETG44aXfjxq0WNAhYiMS35UvjkD+DnuCP4bPsfilyLgr0D06aFV3nRCcsPxTRtwNnB71PxmIDf54aSEu3BnPV7xO5Bk66+nNnpThTfVqPlrvakAG5IXjq9WAuNVdZeIfMLvYPygqjXAlzpYNNubLk9iOL5R1XbgAzhYyysHPglcAHzOx9B8ISKfAWYAx+NOJ/Yrlmh6rsSbRo/IGb7oW5zEWHylqjv9jiEVicjJwDeBp1R1VVfrp6GrcDV8cNet/upjLEknImOAXwGfVNU97hJm/2KnznrOGwmc6Itd4fntSYzFpBgROR14Hler/YzP4fhlIe402heB03HJpl/wanN3A/9U1ejT6/2G1Wh6rtqbRtdciqKWm37GawBwL65Z7yWqutffiPyhqhtwiXauiNTgWmSeqqrv+hxaMtwCTAEme41lwDsI9Z63qWrat8iyRNNz4WszE/HOSUc8j1xu+hER+U/cufjXgStVtV8dcIjIQOBy4BVV3R6xaKE3HZH8qHxxDTAIqOxgWQvuutW9yQzID3bqrIdUdS3uaC36npmrgTWqujn5URk/icincS3vHsHVZPpVkvEEcffPRF/4v8ibfkD/8DlgVtTjWWCr9/ff/QsteaxG0zt+DNwjIvtxX6IrgOuAOUfdyqQdESkHfgtswt20OT3q4u9aVd3jR2zJ5F30/j3wTRGpBxbgmr9/C7hLVftFTb+j9ykie3H30SzwISRfWKLpBap6r4jkAF/FXfBdD9ykqg/7G5nxwSVAPjAGeLOD5TfSf1pd3QpsBj4N/Ah3FP8D4L/9DMokn/UMYIwxJqHsGo0xxpiEskRjjDEmoSzRGGOMSShLNMYYYxLKEo0xxpiEskRjjDEmoSzRGGOMSShLNMYYYxLKEk0CiMg7ItIkIu+JyNhubD9CRPZEj84pIneIyLyI57eKyPYj99D3ichGb9jftCAi5SJSEPE8pd+fiLwuIhv9jiPVpPLn1t3YROQuEfllAkI6yBJNYvwSuB84GdctTbx+Azzkda8eaRqwKOr5QtLTV4D/8juI3iAil+J68R7sdyzGdODHwOdFZEqiXsASTQJ4AxzdDBzADd8aMxE5Czfs7x0dLJ5KP0k0qvqUqr7kdxy95GRggN9BGNMRr4f5vwG/TtRrWKJJEFVtBZYBJ3ij7MXqVuBNVd0SOVNERuKOiBd7z3OACtI00RhjkupB4LxE1Wqs9+YE8ZJLNlAIjMWNWdPVNqOADwP/GTFvI64n4LB3o7qdf1JE3lDVc46y3zNxveae4s2aB/xQVedGvc5LuIOPjwF7gGmquruD/W3EDYewGPg6MAqXVG/B9db7W+BSoAY3qNP3VLXd2zaAG6PjU8AkIAvYCNwD/Dw82qD3Ghsj31dvvw9v/VNxpw7C+3wX+K6qRl4Liyfm6NdeAlzm7WpD9GclIh8FvoMbKG8T8CtV/WNUjNcC38YdWKwDvgn8B5Ab3ldH5dXR/FjfSwflVAi8AhwPXKyqb8dRfqW4o+XzgCG4XpwfAX6kqo0dvV5E7C97+/yOt+1ib/+vRa0bSxwbie+7cT1uWAPBlfsXO1mvy9f21jsZ9/09FTfE+3vAN1X1A295zJ9Nb8cGzAX24b5Xn+1oXz1hNZrEuRmY7v09OcZtLgEyOHxM9a/gupb/O+5LdyOHupqv9v7u9FqGiFyBG+VxNHCb9xgNvOIti3QD7vTcl4E/dfYP6PkI7gt8F64L+ArgcdwPQzvwf3DJ59tejGG3AX8AVuAS6reBRuBnwE3JfB8iciHwBlACfA/4ibfPuV5S607Mh722t88nvWW3cvhnNQv4v8Cj3n6bgD+IyOyIGG/C/Si34JL6G8BjuB/87oi7/EUk23sPU4ArIpJMrOX3CPAhXHncgvscv4k7IOnKhcD/w73n7wHlwAsicnZEfLHGAbF/Nz4BPATU48r9VdzB1ZCo9WJ6be/vucBxuGESfoL7DF+PaDAU02fT27HBwTMwz+MOEHud1WgSQESGA7cDO4ChuETzTAybngHU4cazAdy1Cm+fVwHvqOpfveenAYvDzzuJIxP3T7oNmKmqNd78O3FJ4Pci8pyqtnib5AHXqeq6GGIdAZwYcTQ2EPga8LaqzvHmPYA7SroIN058Fu7I6yFV/UREnHcBu3Cjkt6XjPchIkHgj7ha0dmq2ubN/x3uqPm3wLRuxHzEa4vIUuBK4ClV3Ri17pmqutBb71lczfcq4CkRycD9KK3y1mvy1lvlxRfL5xT5nuMuf6+cHgDOwg1J/WrE/FjKrxy4APiaqv7C2+1d3tH7+BjCHu29bvj/4C/AatyP76mxxhGxv1i+Gxm4a6TzvX22ePMX4moXkWUT62v/AtgLzFDVvd56/wBWAl8Qke8Qw2eToNjClgIfFZFxHTRE6hGr0STG73DV3qu957HWaMbjTnN0dPriRNxpmMjni7vY33RgJPC78I8zgKpWeTGOAGZGrL82xiQDsC6cZDyrvWn46B1VrcP9kwzznrfgjrqiq+aDcKfZCpP4PqbhyvspoFREBonIINwP0d+BqSIyshsxx1OGq8NJxns/m4DduIMTcDWecuDOcJLx3Ansj/E1Dupm+f8RN0z5Z1X1nxHzYyo/XK37AO7H9OpwE29V/ZSqXhBD2KvCScbbbjfwF+BkL4nFGkdYLJ/PdFy53xNx8IL3upHlHtNre3HOAh4MJxnvvazGfW/viOOz6dXYol4rfIA7jl5mNZpe5p32uBL4uqq+IyK7gBNi3LwMqIraXylQhPvwN3lfFnDJ62HveUsn49KHvzAdDZu70puOwZ23BZcUYrUz6nlrJ/to4/ADmmbgchH5CO788jFAqbesswOfRLyPCd70v+l8xMdRuOsJ8cQcTxl2tG4D7toeuKN5iKjhAqhqs4jEVZuJEM97GYMbMRbgdA6v7cRUfqq6VUQ+hztt9hjQJCJv4E6z3n+0azSeFR3MWwMEvPjC341YPkeI7fMZ600PK2NVbRORNRGzYv0OtXrxroleqKqRrUhj+Wx6O7atEc/DB3GDOli3RyzR9CIRKcYdYb8P/MqbvRQ4R0SyVbW5i120c+Q/+yIONQZ4KGrZ/3iPN4BzOtjf0Vq7hV8nMqa2LuKL1NrJ/E6HbPVOl/wVd578LeAd3NH5XNx55s4k4n1keNPv4S7KdmRVN2KOpwzb41g3Wlc/0GHh99md8g/hrjWeDnxGRO4LX58hxvIDUNUHReR5XLP9y3Gn0i7C1XJOjqqtRevofyb82m3xxBGxTVfC3+HcDpZF/n/G+trh1judft5xfDa9HVtH28fzHY6JJZredTuu+vuh8DlRXKK5AHexfGkX2+/k0FFs2MeAT+BaLn3cm3ch8CVcCzXo/DTKRm9aATwdtSz85d9C8pyJ+0e6TVW/fzAQdw2mjKgj9wgbvWlvvo/wPg+o6suH7VBkFjAQV7vobsy9IXyEKh0sG8/hR8htQE7kCl6Mgzh09Bvve9msqneKyFPAFcCdIjLNO2Wz0VvnqOXntVabCixX1buBu73GBT/HXZC/CHcqpzMTOph3jPd+N0S8564+x3iEy+HYqP0FcDWK5d6sjTG+9mZv9sToFxKRO3D/v+8Q22fT27FFKvOm0Wcresyu0fQSETkF+DzwC1WNvHYSTi6xXKfZBAz3LvgB4B1B5gP/UtWXvS9NBrAo/FxV3+9kf+8Dlbgjx+KIWIuBL3jLOts2EcJf5OjTIf+Oe4+dHfgk4n0s8Lb7kvdjGLnPR3AXVlt7EHOk8EFHvP9vS3BJ4vNyePc11wDDo9bd4RZJXsS8Kzj8yLdb70VVdwLfx7WSCvd0EWv5nQC8CXw6Yn/NHLrxuKuj51ne/1Z4/0OAfwNeVdX9ccQRj0W4H+qbRSQ/Yv4cDj+tFNNrq+p23Gd5Q9T3dxwu2Q4h9s+mV2OLeq3wNZvN9DKr0fQCrzXPn3A/Cj+KWhxPonkV+CTunzPywv8s3P0oYTNwX6SjUtUWEfki7ou1wGvBAu68+3DgmvD9LUnyDu488K9FZDTuetS5wPW4U0FFHW2UiPcRtc+F3j4bcf/YY4CPqWqriHQr5ijhJrRf81rHxdICEVVtF5Gbcc3d3xWRu3ENBb7EkaeU/oZrKv28iPwVd/T8WdzBS1hP3sv/w93f8T0ReUhVN8RYfv/CJZr/8l5zKe7awBdxp25e7uC1IjUBz4nIr3FH4LfgEvZXvTKK6XPs4jUOo6ohb59PcajcR+DuMdkXsV48r30r8AIw31uv3SuDKlwrshAxfDYJii3sFFxjiV5PNFaj6R1fxx3t/XsHFzdX4I4cYkk0L+C+gJHt70twPxoLvOcBXIuSLhMNHOwO5yJgO+5msW/jTjmcG9maJxm8I+PLcAn5e8BPcV/6OcDvgeO9I9aOtu319xGxz61ePLfh/tmvUNW/9TTmCA/hflA/ScddCx0txpe8GOtwTXqv8vZTGbXq73HlMg6XcM7BNUpZFrGvnpR/G672mOutG2v5hXDXZv6Iu5fmd7gE+Djus+vquuV7uHtuPourVa0ATlfVg6ehY4kjXqr6LO56UgPulPiVuFrZyqj1YnptdTeYnuut9wPvPb3vvZcd8Xw2vR0bHGwOfSrwXHfKqyuBUKjTa7fGByLyJDBYVc/wOxY/icgmXBPq8/yOJRVJJz0BpJP+8B5ThbibO18Epqrqkq7Wj5fVaFLPL4DTReSIC4f9TDHu/gtjTOLdBLyUiCQDdo0m5ajq2yLyd+AbuPOp/Yq4PpzOxvV2nJAvvTHmEK9RwjW43h8Swmo0qekW4GoR6ahpZ7oL9+X2DK52Z4xJrO/jep+Yn6gXsGs0xhhjEspqNMYYYxLKEo0xxpiEskRjjDEmoSzRGGOMSShLNMYYYxLKEo0xxpiEskRjjDEmoSzRGGOMSaj/D+WG3uIfNodmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as st\n",
    "alpha = 1.0  # Pick them here\n",
    "beta = 1.0   # Pick them here\n",
    "lambda_prior = st.gamma(alpha, scale=1.0 / beta) # Make sure you understand why scale = 1 / beta\n",
    "lambdas = np.linspace(0, lambda_prior.ppf(0.99), 100)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(lambdas, lambda_prior.pdf(lambdas))\n",
    "ax.set_xlabel('$\\lambda$ (# or major earthquakes per decade)')\n",
    "ax.set_ylabel('$p(\\lambda)$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "D. Use the package ``graphviz`` to draw a graphical model representing the situation.\n",
    "Make sure you use the plate notation.\n",
    "Below, we have introduced all the nodes you are going to need, but some nodes should be \"observed\" (observed nodes are filled), a lot of edges are missing, and you need to use the plate notation (see Lecture 5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: hw01_p6_1 Pages: 1 -->\n",
       "<svg width=\"206pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 206.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>hw01_p6_1</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 202,-112 202,4 -4,4\"/>\n",
       "<!-- alpha -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>alpha</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"22.72\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">α</text>\n",
       "</g>\n",
       "<!-- beta -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>beta</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"95.11\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">β</text>\n",
       "</g>\n",
       "<!-- lambda -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>lambda</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"167.5\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">λ</text>\n",
       "</g>\n",
       "<!-- Xn -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Xn</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"162.44\" y=\"-14.8\" font-family=\"Times,serif\" font-size=\"14.00\">X</text>\n",
       "<text text-anchor=\"start\" x=\"172.56\" y=\"-14.8\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\">n</text>\n",
       "</g>\n",
       "<!-- lambda&#45;&gt;Xn -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>lambda&#45;&gt;Xn</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M171,-71.7C171,-63.98 171,-54.71 171,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"174.5,-46.1 171,-36.1 167.5,-46.1 174.5,-46.1\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a1f922fd0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Digraph\n",
    "import os\n",
    "gc = Digraph('hw01_p6_1')\n",
    "gc.node('alpha', label='<&alpha;>')\n",
    "gc.node('beta', label='<&beta;>', style='filled')\n",
    "gc.node('lambda', label='<&lambda;>')\n",
    "gc.node('Xn', label='<X<sub>n</sub>>')\n",
    "gc.edge('lambda', 'Xn')\n",
    "gc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "E. Show that the posterior of $\\lambda$ conditioned on $x_{1:N}$ is also a Gamma, but with updated hyperparameters.\n",
    "Hint: When you write down the posterior of $\\lambda$ you can drop any multiplicative term that does not depend on it as it will be absorbed in the normalization constnat. This will simplify the notation a little bit.\n",
    "<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F. Prior-likelihood pairs that result in a posterior with the same form as the prior as known as conjugate distributions. Conjugate distributions are your only hope for analytical Bayesian inference.\n",
    "As a sanity check, look at the wikipedia page for [conjugate priors](https://en.wikipedia.org/wiki/Conjugate_prior), locate the Poisson-Gamma pair and verify your answer above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "G. Plot the prior and the posterior of $\\lambda$ on the same plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEdCAYAAAA1s6EDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hb1bXw4Z80vdse94oLLBuwcQWDKaaYHkLHJBfSC+EGQirphSSEFJKbLzeBG0JLIPQWejMYMAYb2xi35d57meLpRd8f+8jIssYjaaSRRrPe59FzRqcu7dHMOvucffb2BQIBjDHGmGTxpzoAY4wxmc0SjTHGmKSyRGOMMSapLNEYY4xJKks0xhhjkio71QGkoWZcAq5KdSDGGNNFlAKttJFTfNa8+RCtgUDAF0+x+Hxu2t2L1MrBsXJwrBycTC4Hnw98Pl+ANq6SWY3mUFWBAGV79uyPecOysgIAKivrEh1Tl2Ll4Fg5OFYOTiaXQ3l5MT5f21eB7B6NMcaYpLJEY4wxJqnS9tKZiIwH5gHDVXXzYdYrBm4DLgOKgdnAjaq6qlMCNcYYc1hpWaMREQGeJbpE+DBwBfA94FpgEDBLRMqSF6ExxphopVWNRkSygS8DvwGaolj/ZOB84DxVfdGb9xawDvgqrqZjjDEmhdKtRnMy8FvgD7gaSnvOBqqBV4IzVHUX8CYuARljjEmxdEs0y4ERqvpz3IOT7RkNrFbVlrD5qwFJdHDtWbF+L7srMq/pojHGdERaXTpT1R0xblJG5Cf4q3FPqsbF5/u4zXu0Vm+q4Id3vsuA3kX8+abT8Pt98R6+y8vOzgJiL8NMY+XgWDk4mVwOvnb+3aVbjSZWPiDSc7Y+XHcInaaoIAeAbbtrWLFhX2ce2hhj0lpa1WjiUAmMiDC/xFsWl0Ag9qd3C3P8DB9YyrqtVbwxfxODemXeWUu0MvkJ6FhYOThWDk4ml4PXM0CbunqNRoERIhL+EUd5yzrVtHEDAZivO2ltzcAOjYwxJg5dPdG8DPQAzgrOEJE+wKnAq50dzEljBwBQWdOIbqro7MMbY0xa6lKJRkT6iMhUESkFUNXZwBvAQyLyBRG5BJdgKoC/dXZ8/XoVMmqwe0503oqdnX14Y4xJS10q0QAXAO8CE0PmXQo8A/weuBfYDJypqim5Ix+8fPaB7qSltVPbIxhjTFqy8WgOVdHaGoh7mIDdFXV85bbXAfjWVeM5ZnivRMeX9jL5pmcsrBwcKwcnk8uhvLwYv99XibuVcYiuVqNJe717FDBqkLt89v7yWB8LMsaYzGOJJgmmjO4LwIKVu2husctnxpjuzRJNEkwZ0xcfUFPfzJK1e1MdjjHGpJQlmiToUZzH6GE9AZi7bHuKozHGmNSyRJMkU4/pB8CiVbupa4imf1BjjMlMlmiSZNJRfcnO8tPY3MqClbtSHY4xxqSMJZokKczPZvyocgDmLrPWZ8aY7ssSTRKdcHR/AJat30vl/oYUR2OMMalhiSaJxo0spzAvm0AA3l9uXdIYY7onSzRJlJPtZ/LoPoC1PjPGdF+WaJJsqnf5bN22arbtqUlxNMYY0/ks0STZUUN7UF6aB8CcJVarMcZ0P5Zokszv83HisW6cmjlLttuAaMaYbscSTSeYdqy7fLavuoHlG1MyeoExxqSMJZpO0K9X4YEened8tC3F0RhjTOeyRNNJThrrajUfrNxlXdIYY7oVSzSd5PjRfcnJ9tPY1Mp8tWdqjDHdhyWaTlKYn8OEI3sDMOcja31mjOk+LNF0omljXesz3VTBzorMG87VGGMisUTTiY45ohc9S9wzNW8vtkYBxpjuwRJNJ/L7fQdqNe98tI2WVhvm2RiT+SzRdLJTxrlEs6+6wYZ5NsZ0C5ZoOlmfHgWM8YZ5nv3h1hRHY4wxyWeJJgVOPW4gAIvX7LFxaowxGc8STQpMPKo3RfnZtLQGrKNNY0zGs0STAjnZWZzo9X82+8OtBALW0aYxJnNZokmRU8e5y2c79tWhGytSHI0xxiRPdqoDCCciVwM/AkYA64FbVfX+w6zfB/gtcA6QD8wBblLVVcmPNn6D+xYzcmApa7ZW8caiLYz2GggYY0ymSasajYhcATwAvAxcDLwB3Ccil7exvg94EjgPuBm4BugPzBKRtP/PPX3CIAA+0F1U1jSmOBpjjEmOtEo0wK3AI6p6k6q+pKrXAY8At7Sx/pHANOC7qnq/qj4HXAkMAi7qlIg7YMrovgcaBby92Jo6G2MyU9okGhEZAYwEHg9b9BgwWkSGR9gs35tWh8wLPgVZntgIEy83J+tATwFvLtpqo28aYzJS2iQaYLQ31bD5q72phG+gqouBWcBPRGS0d7/mz8B+4KlkBZpIp413jQJ2V9azZJ31FGCMyTzp1BigzJtWhc0P1lZK29juOuAlYLn3vgG4WFXXxhuIzwdlZQUxb5ednQXEtm1ZWQFjR5bz0Zo9vL1kG6dMHBzzcdNNPOWQiawcHCsHJ5PLwec7/PJ0qtEEQw2/fhScf0gPlCIyBpgL7AIuwbU8+w/wuIickqQ4E+6cE4YBsGDFTnbZ8AHGmAyTTjWaSm8aXnMpCVse6iZveraq7gMQkVeAt4A/ApPjCSQQgMrK2P/hB89UYt32qEGllBXnUrm/kf/MXsNlp42M+djpJN5yyDRWDo6Vg5PJ5VBeXnzYWk061WiC92ZGhc0fFbY81DBgWTDJAKhqAHgbOCbhESZJdpaf6eNdU+c3F22lsaklxREZY0zipE2iUdXVwDog/JmZy4BVqrox0mbAsRGemZmKe9izy5g+fiBZfh/765p4b/mOVIdjjDEJk06XzgB+AdwjIvuAZ3HPwlwJzIQDvQCMxNViqoDbgf8CXhKR3wC1wLXAacFtuoqy4jymjOnL3KU7eG3+Zk4eOwBfe3fYjDGmC0ibGg2Aqt4LfBV3U/8pYDpwrao+7K1yAfAuMNFbfz3ugc3twL3AQ8AQYEbINl3GmZNci7ONO/ezanOkW1LGGNP1+Kzn4ENUtLYGyvbs2R/zhom42XfLffNYt62ayaP78rWLj417P6mUyTc9Y2Hl4Fg5OJlcDuXlxfj9vkqgR6TlaVWjMXDWpCEALNBd7K2qT3E0xhjTcZZo0szk0X0pLcqlNRDgtQWbUx2OMcZ0mCWaNJOT7ecMr1fnNxdupb6xOcURGWNMx1iiSUPTJw4iJ9tPbUMzby/elupwjDGmQyzRpKHSwlymeUM9vzJ/k/XqbIzp0izRpKkZU1yjgF0V9SxctSvF0RhjTPws0aSpAeVFHDfSDanz0vubUhyNMcbEzxJNGjvn+KEArN5SyZot9gCnMaZrskSTxmRoD4b1c51Xv/BepK7ejDEm/VmiSWM+n4/zprpazYKVu9i6uybFERljTOws0aS5ydKXvj1c1xUvWq3GGNMFWaJJc36/j3O9Ws27S7dbtzTGmC7HEk0XMO3Y/pQV5dLSGuDledYCzRjTtVii6QJysrM423uu5s1FW9lf15TiiIwxJnpxDXwmIkfhhkruCwSAXcASVV2VwNhMiOkTBvHsuxuoa2jm1fmbuPiUEakOyRhjohJ1ohGRMbhBya4A+nmzg0NABrx1dgCPAHeq6vIExtntFeRlc8bEQTz37gZenb+Zs6cMpTA/3QZINcaYQ7X7n0pERgK3AZcAdcBbuFEu1wB7cMmmFzAKmAp8Efi6iDwBfE9V1yYn9O7n7ClDeHX+ZmobmnltwWY+cdIRqQ7JGGPaFc0p8TLgI+CzwBOqetiHOUSkCLgcuMHbNr+DMRpPSWEuZ0wcxAvvbeTl9zdy1qTBFORZrcYYk96i+S91pao+He0OvUR0H3CfiHwy7shMROccP5TXPthMTX0zsxZu4fypw1IdkjHGHFa7rc7aSzIi4mtrWSwJykSntCiX6d7AaC++t5GGxpYUR2SMMYeXiOsu+0VkEfABsMCbLlXV1gTs20Rw7glDeX3BFvbXNTFr4RbOPWFoqkMyxpg2JeI5mhnAv4ES4Fu4ZLNfROaKyP+KyOcTcAwTokdxHqeNHwjA83M32HDPxpi01uEajarOAeYE34tIATAeOAH4b1yT6Ls7ehxzsAtOHMZbH7qHN1+dv5kLrQWaMSZNJbRnABHxAycB1wLf9/b/x0Qewzg9ivM4Y9JgwN2rqa233gKMMempwzUaEckCzsI9yPlJ3LM1jwPnqeqCju7ftO28E4Yya+EWahuaeXme9RZgjElPiWgMsBPYDfwLmK6qSxOwTxOFksJcZkwewrNz1vPyvE2cOWkwJYW5qQ7LGGMOkohLZwXASOAq4LsicoOInOQ9uGmS7Nzjh1CYl019Y4uNwmmMSUuJqNEUA2OASd7rCuBXQL6IrAQWqOo10e5MRK4GfgSMANYDt6rq/YdZ34+7H/QFYACwGviVqj4U16fpYgrzczj3hKE8MXstr87fzFmTBtOr1DpjMMakjw7XaFS1VVWXqur9qnqjqp4ClALjgF/jenaOiohcATwAvAxcDLyB62Hg8sNs9ifgx8BfgAuBucCDInJePJ+nK5oxeQhlRbk0t7Ty1NvrUh2OMcYcJCkdZalqAFjuvR6IYdNbgUdU9Sbv/Usi0gu4BXgsfGWvw8/rgS+r6j+82a95wxicC7wQ50foUvJys/jkycO5/yXlnY+2cc6UIQzqU5zqsIwxBkijgc9EZATuXs/jYYseA0aLyPAIm10M1AIHXVpT1dNU9cakBJqmTh43gH69CgkE4PE3rcNsY0z6SJtEA4z2pho2f7U3lQjbjPPWnyEiH4pIs4isEpGrkhVkusrO8nPZqa5586LVu1m5qSLFERljjJNOfcyXedOqsPnV3rQ0wjZ9gKG4ngd+DKzDjYfzkIjsVNVZ8QTi80FZWUHM22VnZwHxbZsIZxw/lFfmb2bV5goen72WW687CZ+vzT5PkybV5ZAurBwcKwcnk8uhvX8zie4Z4JkI856LcvODRuuMMD9SJ525uGTzRVX9u6q+ClwNfAj8LMrjZgyfz8c157uK4apNFbz94dYUR2SMMYmv0bwTYd7bUW5b6U3Day4lYctDVQMtuFZqgGuIICKv4Go2cQkEoLKyLubtgmcq8WybKIN7FTLpqD58sHIX97+wnNGDy8jNyerUGNKhHNKBlYNj5eBkcjmUlxcftlaT0BqNqt4WYd6t0W7uTUeFzR8VtjzUKtxnyAmbn8uhNaNu44rTR5Ll97G3qoGX521KdTjGmG4ubRoDqOpq3D2W8GdmLgNWqWqkx95fxF1auzI4Q0SycU2b30pSqGmvb89CzprsOtx8bu4GKvc3pDgiY0x3lk6NAQB+AdwjIvuAZ4GLcElkJoCI9ME1gV6mqlWq+rqIPA/8WUSKgZXA14DhwKdS8QHSxSdOOoJ3PtrO/romHp+9ls+fPybVIRljuqm0qdEAqOq9uPFrzgGeAqYD16rqw94qFwDvAhNDNrscuAO42dumDzBDVT/onKjTU2F+Dhef4h49emfxNtZtC2/MZ4wxncMXCHTbWxltqWhtDZTt2bM/5g3T7WZfS2srP79nPpt37Wf4gFJ+eO0k/J3Q3DndyiFVrBwcKwcnk8uhvLwYv99XCfSItDypNRoRyRGR15N5DNO2LL+fT884EoB126p4Z/G2FEdkjOmOkn3pzA+cluRjmMOQoT2ZenQ/AB57c42NxGmM6XSJGGHzcDWWtLoH1F1dcfooFq7eTXVtE0+9tY5PzTgq1SEZY7qRRLQ6OxG4HTfSZrgc4JQEHMN0QM+SPC466QgefWMNry3YzLSxAxjWv6T9DY0xJgESkWgWA/NU9anwBSKSD/w2AccwHTRjyhDeWbKdrbtruP+lFfzwmsn4/Z3fD5oxpvtJxKWt+w6zrAn4eQKOYTooO8vPtee4DrDXbatm1sItKY7IGNNdxJxoROTK0Peq+tdItRlvWYuqWqJJE0cN6cHJ4wYA8MTsNeyrth4DjDHJF0+N5kER+UrCIzGd4srTR1FckENdQwsPvbYq1eEYY7qBeBLNvcBfReSHkRaKyIki0m37GUt3xQU5XHm666d03oqdLFy1K8URGWMyXcyJRlW/iLvBf4uI/DE4X5wncMMCTEhciCbRpo3tz5hhPQH450tKbX1ziiMyxmSyuBoDqOr3gW8CN4jIv0TkDuAjXF9kd3JoV/8mjfh8Pj573mhyc/xU7G/k0TdWt7+RMcbEqSOtzv4OvI7rJflLwKPAGFX9mqpuT0RwJnn69CjgslNHAvDmoq0s37AvxREZYzJVPK3OckTkRmANcAawADfIWA5go2x1IWdOGszIQW5A03tfWE59o11CM8YkXjw1mlW4ngB2ARep6mRcreYi4HlvXBjTBfj9Pj573hiys3zsqqjn0TfWpDokY0wGiifRZOEulR2nqs8BqOojuEQzFZglIr0TF6JJpkG9i7jk1BEAzFqwhaXr9qY4ImNMpokn0RypqneramvoTFV9GTgLN7rlnEQEZzrHOVOGMmpQGQB3P7/ceng2xiRUPM2b6w+z7D1cJ5r5HQnKdC6/38cXLhxDbo6ffdUNPPiqPchpjEmchHfjr6rLgWmJ3q9Jrn49C7liumuVPmfJduaviNQZtzHGxC4p48WoqrU+64JOnziIY4b3AuC+F1ewt6rNyqsxxkSt3WECEjQUc0BVz0zAfkwS+X0+vnDBGH7yj/fZX9fEXc8u49szJ9hwAsaYDommRuMHfB182UibXUSP4jw+d/5oAFZsrOCl9zemOCJjTFfXbo1GVad3QhwmjUw4sg+nTxjErIVbeGL2WmRoT0YMLE11WMaYLspqGiaiK88YxcDeRbS0Brjj6SXUWJNnY0ycLNGYiPJysrjuk8eQm+1nd2U9dz+3nEAgkOqwjDFdkCUa06ZBfYr59IyjAFi4ajevfrA5xREZY7qidu/RRENEsoHRwFjvdayqXpSIfZvUOnncAFZsrODdpdt55PXVjBhYysiBZakOyxjThcScaERkGB8nlLHAsYCE7MsH2AX9DOHz+bjmnKNYv72KbXtq+euTS/jpZ6dQWpSb6tCMMV1ENM/RDAW+DUzCJZXQ3pl9wE7gTWBxyGtZvAGJyNXAj4ARwHrgVlW9P8pthwBLgN+p6i/jjcEcLD83m+svGcst989nX3UDdzy9hG/NHE+W3668GmPaF02N5mHgBGA/7p94LjAR2AZcpqpzExWMiFwBPAD8D/AicDFwn4jUqupj7WzrA+4GrB1uEgzsXcTnzx/D355awoqNFTwxe+2BLmuMMeZwojklnYAbTbOnqp7kjT/z30AR8KqIfNP7J58ItwKPqOpNqvqSql4HPALcEsW21+HuE5kkmTK6L+ccPwSAF+ZuZJ71h2aMiUI0ieYD4EVVbQnOUNW/AmOAl4HfA3NFZGxHAhGREcBI4PGwRY8Bo0VkeDvb3oYbJ8ck0eXTRzJ6aA8A/vHsMjZsr05xRMaYdNduolHVaar6ZIT521T1UuBSYBAwX0R+KSLx3iUO1kY0bP5qbyqRNhIRP3Avrib0YpzHNlHK8vu57uJj6V2WT2NzK//vicVU1jSmOixjTBrrcPNmVX1KRF4DfgPcDFwmIl9S1bdj3FWwzWxV2PzgKXNb916+gWs48IkYj9cmnw/Kygpi3i47OwuIb9uupKysgO9/Zgo/vGMOe6sauPOZpfzsiyeQ433+7lIO7bFycKwcnEwuB187N08S0mxIVatV9XrgZKAZeCOO3QRDDX/8PDi/NWw+IiLAL4EvqWplHMc0cTpiQCk3XDkegBUb9nHHk0us5wBjTEQJeWAzSFXnisgE4LtxbB5MFOE1l5Kw5QCISBZwH/Ao8Ir30GiQX0SyVbU5jjgIBKCysi7m7YJnKvFs2xWNHlzGJaeO4MnZa3ljwWZ6FuXwiWnDu105tMXKwbFycDK5HMrLiw9bq2m3RiMiMY0jo6rNqvprb9uzYtnUm4a3mR0VtjxoCK7Z9bW4B0SDL4CfYw+NdooLTxzGtLH9AXjyrXXMXbY9xREZY9JNNDWaF0XkLeB24IXQ1meRiEgOcCHu3smJuOdu2qWqq0VkHXA5ENr44DJglaqGD4yyFZgSYVfzgL/hnqkxSebz+fjMuaPZU1nPio0V3P3ccob0L+Nob6ROY4yJJtFMwCWZZ4DdIvIK8D6wBtiLu4fSCzgSmAqc4b1/CRgfYzy/AO4RkX3As8BFwJXATAAR6YNrAr1MVauA+eE7cLdt2KqqhywzyZGd5ef6S8fy639+wLY9tfzm/nn88isnUVaQ0CuzxpguKprmzUtU9WxgGi55fAL4I/A08BYwG3gK9zzN2biEdIKqnqeqMXVFo6r3Al8FzvH2OR24VlUf9la5AHgX1zOBSSNF+TncdMVxlBXlUlPfzC33vM/eqvpUh2WMSQO+WFsKeTfhJwFHA31wrcR24bqnWaiqh7QO62IqWlsDZXv27I95w0y+2RetjTuq+e2DC6ltaGZAeSHf/69JFBfkpDqslLDvg2Pl4GRyOZSXF+P3+yqBHpGWx5xougFLNB20cXcNv7xnHs0trYwYWMq3rhpPQV73u4xm3wfHysHJ5HJoL9HE9dcvIv1w906OwHW2uQB4RVVr44zTZJCxI3tz41Xjuf3fC1i7tYq/PPER37hi3IEHOo0x3UvMD2yKyCm4bmH+BNwE/BjXSmyDiNyQ2PBMV3XS2AF85lzXq9DyDfv421NLaW7p6ldVjTHxiKdngN97088DQ3G1mpm4Vmh/EpEHExOa6epOPW4gM888EoBFq3fz9/8so6XVko0x3U08l86OBf6gqveFzNsIPCIiXwD+T0TmqOpfEhKh6dLOnjKE+oZmnnp7HfNW7MTv9/HFC8fYoGnGdCPx/LVX4xLLIVT1H8BDuCbKxgDwiWlHcMGJwwB4b9kO/vHcclpbrRGKMd1FPIlmFnB+O8tHxheOyUQ+n49LTx3BeVOHAjB36Q7ues4uoxnTXcSTaP4POElEbmxj+RG47mGMOcDn83H5aSM59/iPk82dzyyzBgLGdAPx3KN5DTcUwO0icglwF24UzmbgNOBG4HsJi9BkDJ/PxxWnjyQry8dz725g/oqdNDe3ct3Fx5KTbfdsjMlU8fx13wI8D2wCTgXux/UKsAK4E9fLcoWIHBvWdb8xBy6jXXyKG5l70erd/PmxD6lvjGtEB2NMF9ChngFEpAeu083x3msCbkjmbFzXNE24xLNYVa/pcLSdw3oG6KBoy+GF9zbw6Kw1AIwYWMo3rjguo7qrse+DY+XgZHI5JKVngCBVrcDd/J8VnCciubgm0MHEM4EEDrNsMsd5JwyjIC+bf76orN1axW8eWMC3rhpPz5K8VIdmjEmghF/aUtVGXJc0CxK9b5N5po8fRHF+Dnc+s5Stu2v41T/nc9MVxzGoT3GqQzPGJIjdgTUpN3l0X75xxXHk5Waxt6qBW/+1AN24L9VhGWMSxBKNSQvHDO/FzZ+aSFlRLrUNzfzh4UU2LLQxGcISjUkbw/qX8MNrJzGgvJDmlgD/98wynn57HTaUhTFdmyUak1Z6lxXwg2smMXqoa7zy9NvruPOZpTQ2taQ4MmNMvCzRmLRTlJ/DN68az6nHDQDg/eU7ue3BheyrbkhxZMaYeFiiMWkpO8vPZ84dzVVnjMIHrNtWxc/vnceqzRWpDs0YEyNLNCZt+Xw+zjl+KDdecRyFedlU1TTy2wcXMmvhFrtvY0wXYonGpL1xI8v58WcnM6h3ES2tAf75knLXs8tpaLT7NsZ0BZZoTJfQr2chP7hmEpOlDwDvLt3OL++fz7Y9NSmOzBjTHks0pssoyMvmuouPZeaZR5Ll97Fldw2/uHc+73y0LdWhGWMOwxKN6VJ8Ph9nTxnC9z41kZ4leTQ0tfCP55bz9/8spa7BeoA2Jh1ZojFd0qjBZfzsc1MYP6o3AO8u3cHP753Hmq2VKY7MGBPOEo3pskoKc/n6ZWO5+qwjyc7ysXNfHbf+cwFPvbXWhok2Jo1YojFdms/nY8bkIfzoWtcqrTUQ4Jl31vPrfy6whgLGpIkODXyWDCJyNfAjYASwHrhVVe8/zPr9caN+ng30wg20dpuqPhpnCDbwWQelqhyamlt4/M21vDxvE+Ae+rzk1OGcM2Uofr+vU2MB+z4EWTk4mVwO7Q18llY1GhG5AngAeBm4GHgDuE9ELm9j/TzgRWAG8BPgUuAD4BEvYZluJCc7i5lnHsl3rp5A77J8mltaeXTWGm791wds2RX7iYMxJjHSqkYjIquB+ao6M2Tew8A4VR0TYf2LgSeB41V1Xsj8F4ABqjo+jjCsRtNB6VAOdQ3NPPbGGmYt3AJAlt/H+VOHceFJw8jJzuqUGNKhHNKBlYOTyeXQZWo0IjICGAk8HrboMWC0iAyPsFkV8H/A/LD5K7x9mW6qIC+ba84RvjNzPH17FNDSGuA/c9bz07vnsXyDDapmTGdK+FDOHTDam2rY/NXeVIB1oQtU9XXg9dB5IpIDXAAsjTcQn+/js49YZHtnyvFsm0nSqRymHjeICUf357HXV/H07LVs31vL7/69kFOOG8hnzh9Dz9L8pB07ncohlawcnEwuB187t0DTpkYDlHnTqrD51d60NMr93AYcCdyaiKBM15eXk8WnzxnNb//7ZGRYTwDe+nArN9z+Js+8tZamZmsKbUwypVONJpgTw28aBecf9r+BiPhwSeYm4Heq+nS8gQQC8V1HzeRrsLFI13LoWZjDd2aO553F23j0jTXsr2vivueX89LcDcw880jGjSxP6PHStRw6m5WDk8nlUF5efNhaTTrVaIKPdIfXXErClh/Ca332IPAdXJL5buLDM5nA7/NxynED+fWXp3LmxMH4fT62763lT49+yO0PL2LzTmudZkyipVOiCd6bGRU2f1TY8oOISCnwCnAl8A1LMiYaxQU5fPrso/jZ56cwxructmTdXn56z/vc/fxyG83TmARKm0SjqqtxN/vDn5m5DFilqhvDtxGRLOBpYCowU1X/J+mBmowyuE8x3545nhsuG8eA8kICAXh78TZuvvNdHnl9NfvrmlIdojFdXjrdowH4BXCPiOwDngUuwtVUZgKISB9cs+VlqloFfBWYDtwJbBKRqSH7Cqjqe0+v+cIAABtRSURBVJ0Yu+mifD4f44/szdiRvZj94TaefnsdVTWNvPj+Rt78cAtnTxnKjMlDKMxPtz8XY7qGtPrLUdV7vfst3wa+CKwFrlXVh71VLgDuAU7H9RpwmTf/K94rVAtp9vlMesvy+zl9wiBOOqY/r36wiefnbqSuoZmn317HK/M2cfbxQzhrkiUcY2KVVj0DpAnrGaCDMqUc9tc18dL7G3l1/mYamtyw0QV52Zw5aTAzJg+mpDD3sNtnSjl0lJWDk8nl0F7PAJZoDmWJpoMyrRyqa91ltNc/2HIg4eTlZHHa+IGcPWUIvdp46DPTyiFeVg5OJpeDJZrYWaLpoEwth+raRl6dv5nXPthMrTeaZ5bfx/Fj+nHuCUMZ0rf4oPUztRxiZeXgZHI5WKKJnSWaDsr0cqhraGbWwi28Mn8TlfsbD8wfM6wnMyYPYdyocvw+X8aXQ7SsHJxMLof2Eo3d1TQmRgV52Zw/dRgzJg9h7tLtvPj+RrbtqWX5hn0s37CPvj0KmD5hEOefPLzd+zjGdAdWozmU1Wg6qLuVQ2sgwLJ1e3l53iaWrNt7YH5utp+TjxvIicf0Y8SAUnzt9TyYobrb96EtmVwOduksdpZoOqg7l8PW3TXMWrCFd5Zso76x5cD8IX2LOfW4gUw9ph9F+TkpjLDzdefvQ6hMLgdLNLGzRNNBVg7uPs6itXt55f2NrN/2cYfk2Vl+Jh7Vm5PHDeDoYb1SMsR0Z7Pvg5PJ5WCJJnaWaDrIysEpKysgEAjwoe7kjYVbeH/FThpCajllxbmceHR/Tjy2P4P7FGXspTX7PjiZXA6WaGJniaaDrByc8HKob2xm/opdvL14Kys3H9wZ+cDeRZwwpi/HH92Pfj0LOz3WZLLvg5PJ5WCJJnaWaDrIysE5XDnsrKhj7tLtzFmynZ37Dl4+tF8xU0b3ZfLovhmRdOz74GRyOViiiZ0lmg6ycnCiKYdAIMD67dW8t2wH7y/fQUXIczkAg/sUMfGoPkw8qg9D+hZ3yctr9n1wMrkcLNHEzhJNB1k5OLGWQ2trgFWbK5i/YhfzV+486GFQgPLSPI4b1Zvxo3ojQ3uSk502o3wcln0fnEwuB0s0sbNE00FWDk5HyqE1EGDtlioWrNzFgpW72Flx8D7ycrIYM6wnY0eWM3Z4L3r3KEhIzMlg3wcnk8vBEk3sLNF0kJWDk6hyCAQCbNldw4erd7No9W7Wbqki/K+2X69Cjj2iF0cf0RMZ2jOthjKw74OTyeVgXdAY08X5fD4G9ylmcJ9iLjjxCKpqG1m6di+L1+5hydo91NQ3s2NvLTv21vLags34fHBE/1LGDOvJ6GE9OHJQD/Jys1L9MUw3ZjWaQ1mNpoOsHJzOKIfW1gAbdlSzZO0elq7by5qtVbS0Hvw3neX3cUT/Eo4c0oOjhvTgyMFlndo7gX0fnEwuB7t0FjtLNB1k5eCkohwaGltYubmCZev3smJjBRt3VBPpT3xg7yJGDSpl5KAyRg4so395If4ktWiz74OTyeVgl86M6UbycrMYO6KcsSPKAaitb2Ll5kpWbqpg1aYK1m+vpqU1wNbdNWzdXcPsD7cBrkfqEQNKOGJAKcO9V4/i3C7ZnNqkH0s0xmSwwvwcxntNogEamlpYv62KVZsrWb2lkrVbq9hf10RdQzNL1+9j6fp9B7YtLcrliP4lDO1XwrB+xQztV0LvsnxLPiZmlmiM6UbycrKQoa5lGrgWbTsr6lizpZJ126pZt62KjTv209zSSlVNI4vX7GHxmj0Hti/Iy2ZwnyKG9C1mcN9iBvcuZlCfIgry7F+JaZt9O4zpxnw+H/16FtKvZyEnHTsAgOaWVrbsqmHDjmrWb69mw/ZqNu/aT1NzK3UNzazaXMmqsL7aykvzGNi7mIG9CxlYXsSA3kUMKC/sdkMimMgs0RhjDpKd5WdY/xKG9S/h1OPcvJbWVrbvrWPTjmo276ph8679bNq5n33VDQDsqWpgT1UDH63dc9C+SgtzGNS3hEF9iuhZnOsltQL69iwgJ9uaXHcXlmiMMe3K8vsZ1LuIQb2LDppfW9/E5l01bPEaF2zdXcPWPTUHus+pqm2iav1elq/fe9B2PqBnaR59exTQx3v17pHvfi4roKQwx+4FZRBLNMaYuBXm53CU93xOqNr6ZrbtrWHb7loqahvZuruGzTuq2bmvjsbmVgLA3qoG9lY1sGJjxSH7zc32U16W716lH796lebRszSfnsV5XaavN2OJxhiTBIX52Ywc6J7RCX1+pDUQoHJ/Izv21rKzoo5dFXXs3OemuyrqqKlvBqCxuZVte2rZtqe2zWOUFubQsySfniV59CjJo0dxLj2K87xXLmXFeZQU5iTt+SATPUs0xphO4/f56FmSR8+SPEYP63nI8tr6ZnZX1rGnsp7dlfXsqfJelfXsrW6gqubjHq2rapuoqm1iw47qwx6vpCiHsqJcSotyKSt005LCXEqLcigtdD+XFOZQUphj942SxBKNMSZtFOZnMzTfPbsTSVNzK/uq69lX3XDgtbe6gYr9DVQEp/sbD3TDE6xBhQ+50Ja83CxKCnIoDnkVBaf52QfeF+ZnU5TvpoV52WRn2WW8w0m7RCMiVwM/AkYA64FbVfX+w6xfDNwGXAYUA7OBG1V1VfKjNcZ0ppxsP317FtL3MCOPBgIBauqbqahuoLKmkcoaN63yXpU1jVTXNlFd66ahfcM1NLbQ0NjC7sr6mOLKy8lyScdLPAV53tR7n5+bRa8ehRTmZRNoaaUgL4v8XDffvbLJzfFnbAOItEo0InIF8ADwP8CLwMXAfSJSq6qPtbHZw8AU4DtANfBTYJaIHKOqlW1sY4zJUD6f70BtZHA76wYCAWobmg8knv21TVTXuZ9r6prZX9fE/romauqbqKl372vrm2huObgDuYamFhqaWg40944vbg4knbycLPJys8j3pnk5H79yc/0f/5yTRV6On9zskJ9zssjJ9pOb7X7OzfaTk51FdpYvZYksrRINcCvwiKre5L1/SUR6AbcAhyQaETkZOB84T1Vf9Oa9BawDvoqr6RhjTEQ+n4+i/ByK8nPo36vtWlKoQCBAY1MrNfVN1NY3H5jWNniv+mbqvJ/rvPf1jc00NLW6efXNNDS1RNgv1DW0UNdw6LJE8AE5OX5ysrxklOUnJ/vg18DeRVx5+qiEXwpMm0QjIiOAkcD3wxY9BlwpIsNVdV3YsrNxtZhXgjNUdZeIvIlLQJZojDEJ5fP5XC0jN4tepdFvF9r6rqW1lYbGFuobW6hrbKG+oZn6phbqG1qob2ymsanlwPuGppYD7xsa3c8NTa3e1L2amltpaGqluaW1zeMHgMamVi9JNkdcZ9n6fZw6biCD+xbHUiTtSptEA4z2pho2f7U3FVxNJXyb1aoafgqwGrgq3kB8vo+/FLHI9lqsxLNtJrFycKwcHCsHpzPKoaU1QGNTy4Ek1NjUSmNzizfPJafGZjdtam6lqdnNb2p26/XrVcjRo3rHfImtvdXTKdGUedOqsPnBtouRzh3KIqwf3CaGcw1jjOn6svw+CrzGCOkknaIJ5sTwYZqC8yPVCX0R1g/Ob7sO2Y5AIL7BiTJ5YKNYWDk4Vg6OlYOTyeVQXl582FpNOjX+DrYQC6+JlIQtD98mUs2lpI31jTHGdLJ0SjTBezOjwuaPClsevs0IEQnPpaPaWN8YY0wnS5tEo6qrcTf7Lw9bdBmwSlU3RtjsZdwY1WcFZ4hIH+BU4NUkhWqMMSYG6XSPBuAXwD0isg94FrgIuBKYCQeSyEhgmapWqepsEXkDeEhEvgvsBX4GVAB/6/zwjTHGhEubGg2Aqt6Le9DyHOApYDpwrao+7K1yAfAuMDFks0uBZ4DfA/cCm4EzVXUfxhhjUs4XCERqtNWttQYCAV88xRJsddHdi9TKwbFycKwcnEwuB58PfD5fgDYqL5ZoDtWMK6xIz+cYY4w5VCnukZKIt2Ms0RhjjEmqtLpHY4wxJvNYojHGGJNUlmiMMcYklSUaY4wxSWWJxhhjTFJZojHGGJNUlmiMMcYklSUaY4wxSWWJxhhjTFJZojHGGJNUlmiMMcYkVbqNR9NlicjVwI+AEcB64FZVvT+lQaWQiIwH5gHDVXVzquPpTCLiB74MfA33fdgBPA38VFWrUxlbZ/JGvr0RVw5DgJXAbar6YEoDSzEReQIYp6rhowlnLKvRJICIXAE8gBvx82LgDeA+EQkfLbRbEBHBDVzXXU9kvgv8BXgO9334A/AZ4NFUBpUC38eNE3UfcCHwCvCAiFyZ0qhSSET+C7gk1XF0Nuu9OQFEZDUwX1Vnhsx7GHfWMiZ1kXUuEcnGncn/BmgCegFDulONxjuL3wP8W1WvD5l/FfAQMEFVF6Uqvs4iIjm4mtwDqvr1kPlvAFmqekqqYksVERkILAFqgAar0ZioicgI3PDSj4ctegwYLSLDOz+qlDkZ+C3uDP57KY4lVUqAfwHhl4dWeNORnRtOyrQApwG3hs1vBPI7P5y0cBfuqsdrqQ6ks3XXSxuJNNqbatj81d5UgHWdF05KLQdGqOpOEflsqoNJBVWtAm6IsOhib7q0E8NJGVVtBT6CA7W8vsDngLOAr6QwtJQQkS8Ck4BjcJcTuxVLNB1X5k3DR+QM3vQt7cRYUkpVd6Q6hnQkIicANwNPqeqK9tbPQJfiavjg7lv9K4WxdDoRGQbcDnxOVXe7W5jdi1066zhvJHDCb3YF57d2YiwmzYjINOBFXK32iykOJ1UW4C6jfR2Yhks23YJXm7sbeF5Vwy+vdxtWo+m4Sm8aXnMpCVtuuhmvAcC9uGa956rqntRGlBqqug6XaGeLSBWuReaJqvpuikPrDNcD44CxXmMZ8E5CvfctqprxLbIs0XRc8N7MKLxr0iHvQ5ebbkREvom7Fv8GcImqdqsTDhHpBVwAvKaqW0MWLfCmgzo/qpS4HOgNbIuwrAl33+rezgwoFezSWQep6mrc2Vr4MzOXAatUdWPnR2VSSUS+gGt59wiuJtOtkozHj3t+JvzG/9ne9CO6h68AU8JezwKbvZ//k7rQOo/VaBLjF8A9IrIP9yW6CLgSmHnYrUzGEZG+wJ+BDbiHNieG3fxdraq7UxFbZ/Juev8VuFlEaoH5uObv3wfuUtVuUdOP9DlFZA/uOZr5KQgpJSzRJICq3isiecC3cTd81wLXqurDqY3MpMC5QCEwDHgrwvJr6D6trm4CNgJfAH6OO4v/KfC7VAZlOp/1DGCMMSap7B6NMcaYpLJEY4wxJqks0RhjjEkqSzTGGGOSyhKNMcaYpLJEY4wxJqks0RhjjEkqSzTGGGOSyhJNEojIHBFpEJG5InJEHNsPEpHd4aNzishtIvJ+yPubRGTroXvo+kRkvTfsb0YQkb4iUhTyPq0/n4i8ISLrUx1Huknn31u8sYnIXSLyhySEdIAlmuT4A3A/cAKuW5pY/Ql4yOtePdQEYGHY+wVkpm8Av0p1EIkgIufhevHuk+pYjIngF8BXRWRcsg5giSYJvAGOrgP244ZvjZqInIob9ve2CIvH000Sjao+paqvpDqOBDkB6JHqIIyJxOth/t/AH5N1DEs0SaKqzcAS4FhvlL1o3QS8paqbQmeKyGDcGfEi730eMJoMTTTGmE71IHBGsmo11ntzknjJJRcoBo7AjVnT3jZDgE8A3wyZtx7XE3DQu2Hdzj8pIm+q6vTD7PcUXK+5U71Z7wM/U9XZYcd5BXfy8WlgNzBBVXdF2N963HAIi4DvAkNwSfV6XG+9fwbOA6pwgzr9WFVbvW19uDE6Pg+MAXKA9cA9wG+Dow16x1gf+rkS/Tm89U/EXToI7vNd4EeqGnovLJaYw4/9IXC+t6t14b8rEfkU8EPcQHkbgNtV9Y6wGK8AfoA7sVgD3Az8N5Af3Fek8oo0P9rPEqGcioHXgGOAc1T1nRjKryfubPkMoB+uF+dHgJ+ran2k44XE/qq3zx962y7y9j8rbN1o4lhPbN+Nq3DDGgiu3L/exnrtHttb7wTc9/dE3BDvc4GbVfUjb3nUv5tExwbMBvbivldfjrSvjrAaTfJcB0z0fh4b5TbnAlkcPKb6N3Bdy/8H96W7ho+7mq/0fm7zXoaIXIQb5XEocIv3Ggq85i0LdTXu8tyNwN/b+gP0fBL3Bb4L1wX8aOBx3D+GVuBbuOTzAy/GoFuAvwHLcAn1B0A98Bvg2s78HCIyA3gTKAN+DPzS2+dsL6nFE/NBx/b2+aS37CYO/l1NAf4f8Ki33wbgbyJycUiM1+L+KTfhkvqbwGO4f/jxiLn8RSTX+wzjgItCkky05fcIcCGuPK7H/R5vxp2QtGcG8L+4z/xjoC/wkoicFhJftHFA9N+NzwIPAbW4cn8dd3LVL2y9qI7t/TwbOBo3TMIvcb/DN0IaDEX1u0l0bHDgCsyLuBPEhLMaTRKIyEDgVmA70B+XaJ6JYtOTgRrceDaAu1fh7fNSYI6q/st7fxKwKPi+jTiycX+kW4DJqlrlzb8TlwT+KiIvqGqTt0kBcKWqroki1kHAcSFnY72A7wDvqOpMb94DuLOks3HjxOfgzrweUtXPhsR5F7ATNyrpfZ3xOUTED9yBqxWdpqot3vy/4M6a/wxMiCPmQ44tIouBS4CnVHV92LqnqOoCb71ncTXfS4GnRCQL909phbdeg7feCi++aH5PoZ855vL3yukB4FTckNSvh8yPpvz6AmcB31HV33u7vcs7ex8RRdhDveMG/w7+CazE/fM9Mdo4QvYXzXcjC3ePdJ63zyZv/gJc7SK0bKI99u+BPcAkVd3jrfccsBz4moj8kCh+N0mKLWgx8CkRGR6hIVKHWI0mOf6Cq/Ze5r2PtkYzAneZI9Lli+Nwl2FC3y9qZ38TgcHAX4L/nAFUtcKLcRAwOWT91VEmGYA1wSTjWelNg2fvqGoN7o9kgPe+CXfWFV417427zFbciZ9jAq68nwJ6ikhvEemN+0f0H2C8iAyOI+ZYynBlMMl4n2cDsAt3cgKuxtMXuDOYZDx3AvuiPMYBcZb/Hbhhyr+sqs+HzI+q/HC17v24f6aXBZt4q+rnVfWsKMJeEUwy3na7gH8CJ3hJLNo4gqL5/UzElfs9IScveMcNLfeoju3FOQV4MJhkvM+yEve9vS2G301CYws7VvAEdzgJZjWaBPMue1wCfFdV54jITuDYKDcvByrC9tcTKMH98jd4XxZwyeth731TG+PSB78wkYbNXe5Nh+Gu24JLCtHaEfa+uY19tHDwCU0jcIGIfBJ3fflIoKe3rK0Tn2R8jpHe9He0PeLjENz9hFhijqUMI61bh7u3B+5sHkJquACq2igiMdVmQsTyWYbhRowFmMbBtZ2oyk9VN4vIV3CXzR4DGkTkTdxl1vsPd4/GsyzCvFWAz4sv+N2I5vcI0f1+jvCmB5WxqraIyKqQWdF+h5q9eFeFL1TV0Fak0fxuEh3b5pD3wZO43hHW7RBLNAkkIqW4M+wPgNu92YuB6SKSq6qN7eyilUP/2BfycWOAh8KW/Y/3ehOYHmF/h2vtFjxOaEwt7cQXqrmN+W0O2epdLvkX7jr528Ac3Nn5bNx15rYk43NkedMf427KRrIijphjKcPWGNYN194/6KDg54yn/AO4e43TgC+KyH3B+zNEWX4AqvqgiLyIa7Z/Ae5S2tm4Ws4JYbW1cJH+ZoLHbokljpBt2hP8DudHWBb69xntsYOtd9r8fcfwu0l0bJG2j+U7HBVLNIl1K676e2Hwmigu0ZyFu1m+uJ3td/DxWWzQp4HP4loufcabNwO4AddCDdq+jLLem44Gng5bFvzyb6LznIL7Q7pFVX9yIBB3D6acsDP3EOu9aSI/R3Cf+1X11YN2KDIF6IWrXcQbcyIEz1AlwrIRHHyG3ALkha7gxdibj89+Y/0sG1X1ThF5CrgIuFNEJniXbNZ76xy2/LzWauOBpap6N3C317jgt7gb8mfjLuW0ZWSEeUd6n3ddyGdu7/cYi2A5HBW2Px+uRrHUm7U+ymNv9GaPCj+QiNyG+/udQ3S/m0THFqrcm4Zfregwu0eTICIyFfgq8HtVDb13Ekwu0dyn2QAM9G74AeCdQRYC76nqq96XJgtYGHyvqh+0sb8PgG24M8fSkFhLga95y9raNhmCX+TwyyFfwn3Gtk58kvE55nvb3eD9Mwzd5yO4G6vNHYg5VPCkI9a/tw9xSeKrcnD3NZcDA8PW3e4WSUHIvIs4+Mw3rs+iqjuAn+BaSQV7uoi2/I4F3gK+ELK/Rj5+8Li9s+cp3t9WcP/9gP8CXlfVfTHEEYuFuH/U14lIYcj8mRx8WSmqY6vqVtzv8uqw7+9wXLLtR/S/m4TGFnas4D2bjSSY1WgSwGvN83fcP4Wfhy2OJdG8DnwO98cZeuN/Cu55lKBJuC/SYalqk4h8HffFmu+1YAF33X0gcHnw+ZZOMgd3HfiPIjIUdz/qdOAq3KWgkkgbJeNzhO1zgbfPetwf9jDg06raLCJxxRwm2IT2O17ruGhaIKKqrSJyHa65+7sicjeuocANHHpJ6d+4ptIvisi/cGfPX8advAR15LP8L+75jh+LyEOqui7K8nsPl2h+5R1zMe7ewNdxl25ejXCsUA3ACyLyR9wZ+PW4hP1tr4yi+j22c4yDqGrA2+dTfFzug3DPmOwNWS+WY98EvATM89Zr9cqgAteKLEAUv5skxRY0FddYIuGJxmo0ifFd3NnelyLc3FyGO3OIJtG8hPsChra/L8P905jvvffhWpS0m2jgQHc4ZwNbcQ+L/QB3yeH00NY8ncE7Mz4fl5B/DPwa96WfCfwVOMY7Y420bcI/R8g+N3vx3IL7Y79IVf/d0ZhDPIT7h/o5InctdLgYX/FirME16b3U28+2sFX/iiuX4biEMx3XKGVJyL46Uv4tuNpjvrdutOUXwN2buQP3LM1fcAnwcdzvrr37lnNxz9x8GVerWgZMU9UDl6GjiSNWqvos7n5SHe6S+CW4WtnysPWiOra6B0xP99b7qfeZPvA+y/ZYfjeJjg0ONIc+EXghnvJqjy8QaPPerUkBEXkS6KOqJ6c6llQSkQ24JtRnpDqWdCRt9ASQSbrDZ0wX4h7ufBkYr6oftrd+rKxGk35+D0wTkUNuHHYzpbjnL4wxyXct8EoykgzYPZq0o6rviMh/gO/hrqd2K+L6cDoN19txUr70xpiPeY0SLsf1/pAUVqNJT9cDl4lIpKadmS7Yl9szuNqdMSa5foLrfWJesg5g92iMMcYkldVojDHGJJUlGmOMMUllicYYY0xSWaIxxhiTVJZojDHGJJUlGmOMMUllicYYY0xSWaIxxhiTVP8fWFm8RpIa6xMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alpha_post = 1.0 # Your expression for alpha posterior here\n",
    "beta_post = 1.0 # Your expression for beta posterior here\n",
    "lambda_post = st.gamma(alpha_post, scale=1.0 / beta_post)\n",
    "lambdas = np.linspace(0, lambda_post.ppf(0.99), 100)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(lambdas, lambda_post.pdf(lambdas))\n",
    "ax.set_xlabel('$\\lambda$ (# or major earthquakes per decade)')\n",
    "ax.set_ylabel('$p(\\lambda|x_{1:N})$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "H. Let's work out the predictive distribution for the number of major earthquakes during the next decade.\n",
    "This is something that we did not do in class, but it will appear again and again in future lectures.\n",
    "Let $X$ be the random variable corresponding to the number of major eathquakes during the next decade.\n",
    "We need to calculate:\n",
    "$$\n",
    "p(x|x_{1:N}) = \\text{our state of knowledge about $X$ after seeing the data}.\n",
    "$$\n",
    "How do we do this?\n",
    "We just use the sum rule:\n",
    "$$\n",
    "p(x|x_{1:N}) = \\int_{0}^\\infty p(x|\\lambda, x_{1:N}) p(\\lambda|x_{1:N})d\\lambda = \\int_{0}^\\infty p(x|\\lambda) p(\\lambda|x_{1:N})d\\lambda,\n",
    "$$\n",
    "where going from the middle step to the rightmost one we used the assumption that the number of earthquakes occuring in each decade is independent.\n",
    "Carry out this integral and show that it will give you the [negative Binomial](https://en.wikipedia.org/wiki/Negative_binomial_distribution) distribution $\\operatorname{NB}(r,\\theta)$, see also the [scipy.stats papge](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.nbinom.html) with parameters\n",
    "$$\n",
    "r = \\alpha + \\sum_{n=1}^N x_n,\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\theta = \\frac{1}{\\beta + N + 1}.\n",
    "$$\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I. Plot the predictive distribution $p(x|x_{1:N})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEdCAYAAAA1s6EDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxcVZn/8U9lgYCBCA3iqCxZ8IkDsqgMIIKArKIxSMDgknEBBBzFMC6giEp0MijihmAEIeHHEmKAoBEJYQmgZJRNEZHHBAiBQRgJIWELhqR+f5xz09W3q7tvVdfpqlS+79erX7frrk+dulXPPfece2+pXC4jIiKSyqBmByAiIu1NiUZERJJSohERkaSUaEREJCklGhERSWpIswNoQa8SEvDKZgciIrKe2BxYSw85paTuzd2sLZfLpXqKpVQKQxWpyiJP5dGVyqNTO5RFqQSlUqlMD2fJVKPpbmW5zIhly16oecERIzYBYMWKlxsd03pHZdGVyqMrlUendiiLjo7hlEo9nwVSG42IiCSlRCMiIkm13KkzMzsWOAMYBSwBprr7pb3MvzXwHeBQYBhwJzDZ3Relj1ZERPrSUjUaMzsauBy4ERgPLABmmNmEHuYvAdcChwOnAR8DXg/camZbDETMIiLSu1ar0UwFZrn75Ph6npltCUwBZleZf0dgH+Dfs1qPmf0VeBgYB8xIH7KIiPSmZWo0ZjYKGA1cnZs0GxhrZiOrLDYsDp+vGPdsHHY0NkIREalHyyQaYGwcem784ji0/ALufj9wK3CmmY2N7TU/Al4A5qQKVEREimulU2cj4jDfFzurrWzew3InAfOAv8bXrwDj3f2RegMplTr7ttdiyJDBQH3LthuVRVcqj65UHp3aoSyyi0570kqJJgs1f31sNn5tfgEzewuhl9li4PPAS8DxwNVmdpi735Eo1m4eeGQZgweHCuKaNSHUnUfp7J2ISCslmhVxmK+5bJabXinrNHCIuy8HMLP5wB3A94F31BNIuVz7Vbr3PPgUw4YNBWDVqtUAbNuxaT2bbwvtcLVzI6k8ulJ5dGqHsoh3BuhRK7XRZG0zY3Ljx+SmV9oeeDBLMgDuXgZ+C+zU8AhFRKRmLZNo3H0x8CiQv2bmKGCRuy+tthiwc5VrZvYiXOwpIiJN1kqnzgDOAi4xs+XAXMK1MMcAE2HdXQBGE2oxK4FzgY8Srrf5b0IbzSTg3dkyIiLSXC1TowFw9+nAiYTbycwB9gcmuftVcZYjgIXA2+L8SwgXbD4FTAdmAtsCB1csIyIiTdRqNRrcfRowrYdp0wkJpXLcXwk1HxERaUEtVaMREZH2o0QjIiJJKdGIiEhSSjQiIpKUEo2IiCSlRCMiIkkp0YiISFJKNCIikpQSjYiIJKVEIyIiSSnRiIhIUko0IiKSlBKNiIgkpUQjIiJJKdGIiEhSSjQiIpKUEo2IiCSlRCMiIkkp0YiISFJKNCIikpQSjYiIJKVEIyIiSSnRiIhIUko0IiKSlBKNiIgkpUQjIiJJKdGIiEhSSjQiIpKUEo2IiCSlRCMiIkkp0YiISFJKNCIikpQSjYiIJKVEIyIiSSnRiIhIUko0IiKSlBKNiIgkNaTZAeSZ2bHAGcAoYAkw1d0v7WX+QcDpwKeAfwEWA99295npoxURkb60VI3GzI4GLgduBMYDC4AZZjahl8V+AHwNOA94H/A/wBVmdnjaaEVEpIhWq9FMBWa5++T4ep6ZbQlMAWbnZzaz0cBngBPc/edx9M1m9mbgMOA3AxCziIj0omVqNGY2ChgNXJ2bNBsYa2Yjqyw2HngJ6HJqzd3f7e6nJAlURERq0ko1mrFx6Lnxi+PQgEdz03aJ8x9sZlOBneI8Z7j7VfUGUirBiBGb1LTMsGFDGTSotO5/qH0d7WTIkMHAhl0GlVQeXak8OrVDWZRKvU9vmRoNMCIOV+bGPx+Hm1dZZmtgO+BiQhvNYcA9wEwzOyBFkCIiUptWqtFkObHcw/i1VZbZiJBs3u/ucwHM7GZC7egbwK31BFIuw4oVL9e0zKpVq9fVZFatWg3Uvo52kh2dbchlUEnl0ZXKo1M7lEVHx/BeazWtVKNZEYf5mstmuemVngfWEHqpAeDuZWA+4bSaiIg0WSslmqxtZkxu/Jjc9EqLCO9haG78RnSvGYmISBO0TKJx98WEhvz8NTNHAYvcfWmVxW4gnFo7JhthZkMIbTV3JApVRERq0EptNABnAZeY2XJgLjCOkEQmApjZ1oQu0A+6+0p3v8XMrgd+ZGbDgb8BJwMjgQ834w2IiEhXLVOjAXD36cCJwKHAHGB/YFJFV+UjgIXA2yoWmwD8FDgtLrM1cLC73zMwUYuISG9K5bKaMnKeW7u2PGLZshdqWmjOHY9063U2ft9RDQ9ufdEOPWkaSeXRlcqjUzuURUfHcAYNKq0AXlttekvVaEREpP0o0YiISFJKNCIikpQSjYiIJKVEIyIiSSnRiIhIUko0IiKSlBKNiIgkpUQjIiJJKdGIiEhSSjQiIpKUEo2IiCSlRCMiIkkp0YiISFJKNCIikpQSjYiIJKVEIyIiSSnRiIhIUko0IiKSlBKNiIgkpUQjIiJJKdGIiEhSSjQiIpKUEo2IiCSlRCMiIkkNqWchM3szsBPwOqAM/AN4wN0XNTA2ERFpA4UTjZm9BTgROBrYJo4uxWE5zvM0MAuY5u5/bWCcIiKynuoz0ZjZaOBs4EjgZeAOYCHwMLCMkGy2BMYAewHHAZ81s2uAL7v7I2lCFxGR9UGRGs2DwJ+BjwPXuPuLvc1sZq8BJgCfi8sO62eMIiKyHiuSaI5x9+uKrjAmohnADDP7QN2RiYhIW+iz11ktSaaRy4qISHuouXuzmR2TIhAREWlP9VxHc4WZfbrhkYiISFuqJ9FMB843s69Wm2hme5vZHf2KSkRE2kbNicbdjwO+A0wxs+9n4y24BvgtsHvjQhQRkfVZXXcGcPfT48WZ3zOzrYEXgE8SLtycBpzVuBBFRGR9VleiiS4EjgA+TEgwM4Gv9fcCTTM7FjgDGAUsAaa6+6UFl90WeAD4rrt/qz9xiIhIY9TT62yomZ1CuDPAgcC9hEQzFHi8P8GY2dHA5cCNwHhgAeF6nAkFli0BFwOb9ycGERFprHo6AywCziXcSHOcu7+DUKsZB1xvZsP7Ec9UYJa7T3b3ee5+EuHeaVMKLHsSMLYf2xYRkQTqSTSDgeOBXd391wDuPouQaPYCbjWzrWpdqZmNAkYDV+cmzQbGmtnIPpY9O8YlIiItpJ5Es6O7X+zuaytHuvuNwEHASODOOtab1UY8N35xHFq1hcxsEKHL9Sx3v6GO7YqISEI1dwZw91W9TPu9me0LzKsjlhFxuDI3/vk47Knt5fOEjgPvr2ObVZVKMGLEJjUtM2zYUAYNKq37H2pfRzsZMmQwsGGXQSWVR1cqj07tUBalUu/TG/6Ezfgcmn3qWLTLs22qjF+bG4+ZGfAt4Hh3X1HHNkVEJLEiz6N5j7vfXMtK3f3xuOxB7n5TwcWyRJGvuWyWm57FNZhwl+hfAPPNrPK9DDKzIe7+ai1xZ8plWLHi5ZqWWbVq9bqazKpVq0PANa6jnWRHZxtyGVRSeXSl8ujUDmXR0TG811pNkVNnN8RbypwL/Mbd1/Q2s5kNBd5HOKW1N7BRwViztpkxhOffUPG6cnpmW2DP+DcpN+2b8a+PCp2IiKRWJNHsTkgyvwSeMbP5wB8I19E8S+cTNnck9Do7ML6eB+xWNBB3X2xmjxIemnZtxaSjgEXuvjS3yJPAHlVWdRdwAeGaGhERabI+E427PwAcYmZ7AycDHwCOpXpbykrgGuACd7+rjnjOAi4xs+XAXEKX6WOAiQDxdjejgQfdfSVwd34FodmGJ9292zQRERl4hXuduftCYGFsG3k78K/A1oSE8w/CrV/uy3d7roW7TzezjYEvAMcBjwCT3P2qOMsRwCXAAYS7BoiISIsrlcv5iskG77m1a8sjli17oaaF5tzxSLfOAOP3HdXw4NYX7dDA2Ugqj65UHp3aoSw6OoYzaFBpBfDaatPruqmmmW1DOKW1A+HOzfcC8939pTrjFBGRNlVzookXZF4PbErXXl3LzGyKu/+oUcGJiMj6r54LNs+Jw08C2xFqNRMJvdB+YGZXNCY0ERFpB/WcOtsZ+J67z6gYtxSYZWafAn5mZne6+3kNiVBERNZr9dRonicklm7c/eeEB6Cd2J+gRESkfdSTaG4F3tvH9NH1hSMiIu2mnkTzM+Cd8Smb1exAuGpfRESkrjaam4FXgXPN7EjgIuCeOO7dwCnAlxsWoYiIrNfqSTRTgF0J9zHbL/5VXvV5H/Ccme0MPFTvHZRFRKQ91PPgs69n/5vZawk33dwt/u0O7AJcRkg+q83Mgfvd/WMNiVhERNYrdd0ZIOPuzxEa/2/NxpnZRoQu0Fni2Z0GPv1SRETWL/1KNNW4+z8Jt6S5t9HrFhGR9U/DH+UsIiJSSYlGRESSUqIREZGklGhERCQpJRoREUlKiUZERJJSohERkaSUaEREJCklGhERSUqJRkREklKiERGRpJRoREQkKSUaERFJSolGRESSUqIREZGklGhERCQpJRoREUlKiUZERJJSohERkaSUaEREJCklGhERSUqJRkREklKiERGRpJRoREQkqSHNDiDPzI4FzgBGAUuAqe5+aS/zvx6YAhwCbAk4cLa7/yJ9tCIi0peWqtGY2dHA5cCNwHhgATDDzCb0MP/GwA3AwcCZwAeBe4BZMWGJiEiTtVqNZiowy90nx9fzzGxLQo1ldpX5Dwd2Bf7N3e+K4+ab2XbAl4ErUwcsIiK9a5kajZmNAkYDV+cmzQbGmtnIKoutBH4G3J0b/1Bcl4iINFkr1WjGxqHnxi+OQwMerZzg7rcAt1SOM7OhwBHAXxLEKCIiNWqlRDMiDlfmxj8fh5sXXM/ZwI6ENp66lEowYsQmNS0zbNhQBg0qrfsfal9HOxkyZDCwYZdBJZVHVyqPTu1QFqVS79NbKdFkoZZ7GL+2t4XNrERIMpOB77r7dY0NT0RE6tFKiWZFHOZrLpvlpncTe59NByYSksyX+hNIuQwrVrxc0zKrVq1eV5NZtWo1UPs62kl2dLYhl0EllUdXKo9O7VAWHR3De63VtExnADrbZsbkxo/JTe/CzDYH5gPHAJ/vb5IREZHGaplE4+6LCY39+WtmjgIWufvS/DJmNhi4DtgLmOjuP0weqIiI1KSVTp0BnAVcYmbLgbnAOEJNZSKAmW1N6Lb8oLuvBE4E9gemAY+b2V4V6yq7++8HMHYREamipRKNu0+P7S1fAI4DHgEmuftVcZYjgEuAAwh3DTgqjv90/Ku0hhZ7fyIiG6KW+yF292mEGkq1adMJjf7Z6wMHJioREalXy7TRiIhIe1KiERGRpJRoREQkqZZro5H+eeix5d3Gjd1+iyZEIiISKNG0mYeWKtGISGvRqTMREUlKiUZERJJSohERkaSUaEREJCklGhERSUqJRkREklKiERGRpJRoREQkKSUaERFJSolGRESSUqIREZGklGhERCQpJRoREUlKiUZERJJSohERkaSUaEREJCklGhERSUqJRkREklKiERGRpJRoREQkKSUaERFJSolGRESSUqIREZGklGhERCQpJRoREUlKiUZERJJSohERkaSUaEREJCklGhERSUqJRkREkhrS7ACk/Tz02HJeM/wlAF584RUAxm6/RTNDEpEmUqKRhnto6XKGDRsKwKpVqwElGpENWcslGjM7FjgDGAUsAaa6+6W9zD8cOBs4ChgO3A6c4u6L0kcrIiJ9aak2GjM7GrgcuBEYDywAZpjZhF4Wuwo4GvgyMAl4I3CrmY1IG62IiBTRajWaqcAsd58cX88zsy2BKcDs/Mxm9i7gvcDh7n5DHHcH8ChwIqGmIyIiTdQyNRozGwWMBq7OTZoNjDWzkVUWOwR4HpifjXD3fwC3ERKQiIg0WSvVaMbGoefGL45DI9RU8sssdvc1VZb5UL2BlEowYsQmNS0zbNhQBg0qrfsfal9HI2TbrjTQcbRKWTzwyLJu43Ye1THgcQAMGTIYaE45tCKVR6d2KItSqffprZRosjaVlbnxz8fh5j0sk58/W6ba/MnsNKqDwYNDBXHNmrUDuelucTRbq5SFiLSGVko0WU4s9zC+2i9Wqcr82fi6f+HKZVix4uWaltm2Y9N1RyTZsrWuoxG27di027iBjkNl0V2+PDZ0Ko9O7VAWHR3De63VtEwbDbAiDvM1kc1y0/PLVKu5bNbD/CIiMsBaKdFkbTNjcuPH5KbnlxllZvlcOqaH+UVEZIC1TKJx98WExv78NTNHAYvcfWmVxW4EXgsclI0ws62B/YCbEoUqIiI1aKU2GoCzgEvMbDkwFxgHHANMhHVJZDTwoLuvdPfbzWwBMNPMvgQ8C3wDeA64YODDFxGRvJap0QC4+3TChZaHAnOA/YFJ7n5VnOUIYCHwtorFPgj8EjgHmA48AbzH3ZcPSNAiItKrUrlcrdPWBu25tWvLI5Yte6HmBduh90ijqCy6Unl0pfLo1A5l0dExnEGDSisITRndtFSNRkRE2o9qNN2tLZfLpXqKJetHriJVWeSpPLpSeXRqh7IolaBUKpXpofKiRNPdq4TCqnbHARER6W5zwkXyVTuYKdGIiEhSaqMREZGklGhERCQpJRoREUlKiUZERJJSohERkaSUaEREJCklGhERSUqJRkREklKiERGRpJRoREQkKSUaERFJqtWesLneMrNjgTOAUcASYKq7X9rUoJrAzAYBJwAnE8riaeA64Ovu/nwzY2sFZnYNsIu7j2l2LM1iZvsB/0V4gOFzwNXA6e5e+0Og1nNmdiJwCrAd8DBwtrtf3tyoGk81mgYws6OBy4EbgfHAAmCGmU1oZlxN8iXgPODXhLL4HvDvwC+aGVQrMLOPAkc2O45mMrO9gPnAU4RHtZ8FfBS4qJlxNYOZnUB45PyvgQ8ANwGXxd+TtqK7NzeAmS0G7nb3iRXjriIcub6leZENLDMrAcuAK939MxXjPwTMBHZ39z82K75mMrM3AA8ALwKvbKg1GjO7Lf67v7uX47jPAKcCb3X3l5oW3AAzszuBVe5+YMW424E17n5A8yJrPNVo+snMRgGjCdX/SrOBsWY2cuCjaprNgMuAK3LjH4rD0QMbTku5iFDjvbnZgTSLmW0F7AtckCUZAHf/ibuP3pCSTDQMyJ9OXgZ0NCGWpNRG039j49Bz4xfHoQGPDlw4zePuK4HPVZk0Pg7/MoDhtAwzOw54O7ATcE6Tw2mmtwIl4NlY438f4UGDVwCnuvvLzQyuCX4IXBhPlc0DDiGUyVeaGlUCSjT9NyIO80/kzI5UNh/AWFqOme0JnAbMcfeH+pq/3ZjZ9sC5wCfc/Rkza3ZIzbR1HE4HrgXeD+wKfAvYBPh4U6JqniuBA4FZFeNmuPt3mxRPMko0/Ref+E2+sSsbv3YAY2kpZrYPMJdQozuuyeEMuNhmdTFwvbvnT61uiDaKwzsr2vBuieV0jpmd5e6PNCm2Zvgl8E5C+9S9wJ7AmWa20t2rnRlYbynR9N+KOMzXXDbLTd+gxA4A04G/AYe5+7LmRtQUnwF2Ad5qZtl3rQQQX6+pbKvYAGS1/Otz4+cReie+FdggEo2ZvRM4lFDTnR5H32ZmzwHTzOxCd/9z0wJsMHUG6L+sbSbfi2hMbvoGw8xOJZwWWAjs5+5/b3JIzTIB2Ar4O7A6/k0idIpYTej2vSFZFIcb58ZnNZ0NKeluH4e/y42/PQ7/dQBjSU6Jpp/cfTHh1FD+mpmjgEXuvnTgo2oeM/sU4eh0FqEms0HW6KJPA3vk/uYCT8T/f9W80Jrir8BjwMTc+KxTwMIBj6h5sgPQ/XLj947DJQMXSnq6jqYBzOzjwCXATwg/JOOAk4CJ7n5VE0MbUGb2OkLS/QfhIrxXc7MsdvdnBjywFmJm04F3bcDX0XyIUNu9gnBq9e2EizbPc/f/bGJoA87MrgUOAs4E7gPeEf//rbu/t5mxNZraaBrA3aeb2cbAFwiN3o8AkzakJBMdBmxKOC1wR5XpHyNcZyMbKHe/ysxeIfygzgX+j5BopjY1sOaYCHwdmAy8jlCLOQc4u4kxJaEajYiIJKU2GhERSUqJRkREklKiERGRpJRoREQkKSUaERFJSolGRESSUqIREZGklGhERCQpJZomMrM3mtkz+adwmtnZZvaHiteTzezJBm2zFNe/zMxeNLOTG7Hefsb0OjN7TcXrJWa2oIkh9crMFpjZkmbHUU2jy66V32uj5ffDGpdt2X223tjM7CIz+14jYlCiaa4fADPdPf8Ezt0J9z6qfH1vg7Z5BPAlwg0MT6HJjxY2s8MJNxjcuq95pSm+DXy+2UGkpv2wqrOAE81sl/6uSPc6axIz24/wiONRVSbvBlxT8Xp3whMJGyHbaU5vkedd7Am8ttlBSHXuPr/ZMQwQ7Yc57r7UzK4Evg+8pz/rUo2meSYDd7j745UjzexNhKOqP8bXGwNjaVyNJnv2x/O9ziUiEu6yfWB/azWq0dTBzDYhPMRpLbCju79SMe0i4BPAR9x9Zg/Lb0t4XvqpFeOW0PkwJICFuefLX2tmt7n7/r3EtS/hbrB7xVF/AL7h7rdX2cajZvaYu+/Qy/r2JlSfs/UtBM5w98r2oxLhuSufBN4CDCXchfYS4DvZEyTjtucTDm4+AjwD/AnIbof+aP79mdmHga8SHiL3GHCuu/80F+PRwFcIyfhh4DTgP4Bh2britpfkyy4/vuh7qVJOwwmnIHcCDnX339VQflsQjhgPBLYhPKtmFvBNd19VbXsVy34IOB2w+N4/W2Weou99Cd0/n92BXwA7ZPtJnO8G4Ldx26OBx4EfuPtPcts4HPgGsDPwNHAuobZ+UB/7XS3b6LWMzeww4DfAL9z9mIrlfgYcDxxOuIty9hC6bvthlfj6LPcisVXMtyfhe7s34Tflf4DTsjMOteyXjY6N8CC2ZwnfqRN6KJI+qUZTB3d/mbBjbAusa0w3s6nAp4DP9pRkosOAwcCvK8Z9nnAb/V8RdqKP0Xlb/RXx/2/3tEIzGwcsALYDpsS/7YCb47RsG9kpuMn0cu7dzA4GbgNGAF8DvhXXd3tMaJkpwAXAg4TE+RVgFfDfhKdJVjqW8ENzCnBhXGdlPJXvbw/gx4QfulOBV4ALzGx8RYyTCD/KqwntTrcBswk/+PWo5b1kMWwU38MuwLiKJFO0/GYRHvx1IeHRzwsIyfJHvQUan4E0E3iJ8N5vIdx2f5ua3nFXXT4fd/9HD/MdHuObTfjcXgTOM7N1z1Axs/cR9uWNCOU4m3AL/CMLxlJkG32WsbvfAMwAjjazQyuWOx74aZw+jZ73wy6KlnvRzz/+fzvhiZrfjfPtBCwwsx3ibIX2y0bHFsvvVULSP7ynMilCNZr6TSfslKeb2YWE59CcBnzd3c/vY9l3Eb44656P7u5zAMzsg8Cd7n5ZfP1O4I/Z62ri8+d/Avwv8A53XxnHTwMeAM43s9+4+xwz243wZZ/j7kt6WN8g4KeEGtG73X1NHH8e4ZTej4DdzWwo4Yhpprt/vGL5iwjPGTmK8CXPbAIc4+4PV8x7fw/xbALs6+73xvnmEh6q9kFgjpkNJnwxH4rzvRLneyjG9zA1qOO9ZOV0OeEpiUe6+y0V44uU3+sID776orufE1d7UTyCrdZ2l213MOGZJXfF9a+O4+8lHOXWq9vn04Ntgd3c/f643WuBJwk1oevjPD8g7N/vjAdmmNnvgDkUO23b6zaKlnFc12TgEEKi2gu4iLB/fAHA3Rf2sh+uU7Tca4ztHGAZ8HZ3Xxbn+zXhaaQnm9lXKbBfJootcz/wYTMbWaXjUiGq0dQpfkCnEdpT5hBOC/zY3c8qsPgowmmLaqdidiWcUqp8/cc+1vc24E2EpxSurIjxOeA84I2Ep/cVtXuMcQ6whZltZWZbEX6IfgXsZmZvijvzNnSvUm8FrASG58YvLvAjlvlblmTie3mM8OTO18dRexAeFjWt8tQl4eh0ecFtrFPHe4HwhZ0AnODu11eML1R+hJrqC4QflKMsdq1190+6+0G9hPs2wnu/JPtBif4fdbz3CkU/H88SQHzxFOHU2OsB4vn80YQaw8sV811H+AEtotdtULyMcfflwImEU7B/IHxXJrn7iwVjyRQt90KxxQONPYArsiQT4/0b4ft6dg37ZUNjy20rOyAeSZ1Uo+kHd58bjxjeQ6iynlJw0Q7gucoR8Vz9ZoQP87H44QO8Fbgqvl7t7iuqrC/bAbzKtOyLvT3Fn8k+Og6/G/+q2ZbQnvBP4Agz+wDhvPCOwBZxnvyBzP8V3H5P875MZ2eG7eLwkcoZ3P2fZlZTbaZCLe9le0ItFmAfutZ2CpWfuz9hZp8mnDabDbxiZrcBVwOX9tJGs0Mcdnmf7r7GzBb1sEwRRT+faqfUXiGcDoZQbhDaMfOc7kfM9Wyjln0Ud/+lmV1DqBFf4O53Foghb4c47Kvci8b2KlCiSjm5e+XlDUX2y0bH9kTF6+zgdasq8xaiRNMPZnYM4Zw2wPM9NRZXsZbuP1z30dlQn2/f+WH8uw3Yv8r6Sr1sK9vOPwvGBp1f5q8RGiareSie4rmMcG7/t8CdhBrF7YTzw3lraohhbQ3z5vXaiF4he5/U8V7KwEmEJHOcmc3I2mcoWH4A7n6Fmd1A6Op+BOFU2iGEWs6eudpa5bYBhlWZVvQsxeAq44p+Pn19NkPjsFrsRT+bvrZRuIwBYm0xS3CHmdlr6qjRFC33orFlvX16fK817JeNjq3a8rV8f7tQoqmTmR1CqJZeS2iM/qSZfd/di5waeJrOI/LMR4CPE3phZT1gDgY+R+ihBj2fFlkSh2OB6/KhxuHjFJet7wV3v6nLysz2ALYk1C72JXwBprj7mRXzDCHU2rrUNhosO0qzKtNG0fUocQ2wceUMMcat6DwCrPW9LHX3aWY2BxgHTDOz3eNpiyVxnl7Lz0Jvtd2Av7j7xcDFsXPBdwi140MIpzPysljenFt3iXBk+5ca33ujVcZ3Y27ajjTGkjjsax/NTCWUzRcJ5f6/sH4AAAReSURBVDuV8N2qRdFyLxrb0jh6TH5DZnY24ft+J8X2y0bHVqkjDp/Ox1mU2mjqYKE74jXA7wgJ4gzCUcnUgqt4DHhDbMADIB4Nbwr83t1vijvBYOC+7LW739PD+u4B/k44Ct68Is7NCb3i/h7nKeruuMzn4o9h5fpmERoXX6VzB3wwt/zx8b0UOZDJjpJq3Rf/RPihPNG63r5mAvCG3LxPhUm2ScW4cXQ9+qvrvbj708CZhJ5CX4iji5bfzsAdhJ6K2fr+SeddIXo6gryP8INxkpltWjF+It1PbxR57412N+HA5lMWrgODEMRehLaERm2jSBljZvsQevT9LHa6uBj4DzN7V8X6iuyHRcu9UGzu/iRhPz42970dSTjQ2Ibi+2VDY8ttK2uzWUqdVKOpkZm9hdAt+W/A+Hhq42Ez+znhR2+filMoPbmFcK3NznRt+N+D0Jst83bCjtErd19tZp8l7Ch3xx4pENoQ3gBMcPfCp6Jy67s3rm8VYefennCN0Ktmdifh/O33zWw7QrvTAcCH4vybFdhcdi7+ixZ6xv2yYIxrzewkwmex0MwuJjQUf47upwmvJHSVvsHMLiMcQZ5ASPiZ/ryXnxCucfiamc1090cLlt/vCYnm23Gb9xPOj3+WcPripirbwt3Lcf1zKt77GwnXOjxbx3tvqPjZnEp4/3ea2aWETjOnEE6nFT3F3Ns2iu6jw4CfE/az0+LiXyacqrzYzHaNHRb63A+LlnvR2OLsk4F5wF1xvrWEz/85Qi+yMgX2y0SxZfYidBSpO9GoRlOD+EHfSOgtdHhlDy/CxU8vE6rlfZlH2KEq+9OPIPwI3B1flwjnlPtMNADufjXhVMuThGt8vkLoDnxA1nW6FhXre4JwPncKYYcf5+5XxnmeJpzqezjO81+EnXUicD6wk5n1dV3HTMIP6icIX6xaYpwfY3yRcE3BB+N6/p6b9XxCmYwk/OjuT+jK+kDFuup+L7EH4smEWsL5cVyR8isTfvB+SriW5jxCEria8Ln12K7m7nMJbTovE2rSRxJqRvlTt32+9xTcfTbhx3AI4TvxYcI1IHdTve2mnm30WcaEC0YN+M/YC5PYw+tLhNN4WS/RQvth0XIvGBvufishcTxB+JxOI5x92Mfdn6plv2x0bLCuO/TehIte61Yql/t9cCF1sHBdwNbu/q4+Z5aaWA9Xw8vAiKeEt/QqF3ya2Z+B5e6+38BHJrWycHHnjYRrmv7U1/w9UY2mec4B9jGzbg2BIuu5wcD/mln+dkE7E9qy8rc5kdY1CZjfnyQDaqNpGnf/nZn9inC++PhmxyPSKPFaplmEbt9lwqmgfyGcXnwGaMgzTiSt2ClhAuHOF/2iGk1zfQY4ysxG9zmnyPrleOCbhJuF/piwr98M/Ju759vQpDWdSbjzxl39XZHaaEREJCnVaEREJCklGhERSUqJRkREklKiERGRpJRoREQkKSUaERFJSolGRESSUqIREZGk/j8kWi1ehJhfsQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "r = 1.0 # Your expression for r here\n",
    "theta = 0.2 # Your expression for theta here\n",
    "X = st.nbinom(r, 1.0 - theta) # Please pay attention to the fact that the wiki and scipy.stats\n",
    "                              # use slightly different definitions\n",
    "fig, ax = plt.subplots()\n",
    "xs = range(0, 10)\n",
    "ax.vlines(xs, 0, X.pmf(xs), colors='b', lw=5, alpha=0.5)\n",
    "ax.set_xlabel('$x$ (# of earthquakes during next decade)')\n",
    "ax.set_ylabel('$p(x)$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J. What is the probability that at least one major earthquake will occur during the next decade?<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "K. What is the probability that at least one major earthquake will occur during the next two decades?<br>\n",
    "**Answer:**\n",
    "<br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L. Find a 95\\% prediction interval for $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here and print() your answer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "M. Find the $\\lambda$ that minimizes the absolute loss (see lecture), call it $\\lambda^*_N$.\n",
    "Then, plot the fully Bayesian predictive $p(x|x_{1:N})$ to $p(x|\\lambda^*_N)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here and print() your answer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-End-"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
